|
Oasis3 4.0.2
|
|
Oasis3 4.0.2
|
00001 SUBROUTINE filling (pfild, pworka, pworkb, pworkc, 00002 $ plon, plat, klon, klat, kmask, kmesh, 00003 $ kunit, cdfic, cdmet) 00004 C**** 00005 C ***************************** 00006 C * OASIS ROUTINE - LEVEL 3 * 00007 C * ------------- ------- * 00008 C ***************************** 00009 C 00010 C**** *filling* - Fill up (and smooth or not) regional field with climatology 00011 C 00012 C 00013 C Purpose: 00014 C ------- 00015 C Performs smooth or abrupt blending of a regional data set with a 00016 C global one for a Sea Surface Temperature or Sea Ice Extent field. 00017 C The frequency of the global set can be interannual monthly, 00018 C climatological monthly or yearly. Sea-ice transitions are treated 00019 C specifically: the ice is supposed to form or disappear on the first 00020 C day of the month. 00021 C 00022 C N.B : The SST climatology must be in celsius. If not 00023 C one has to change the value of local variable zmerg 00024 C 00025 C** Interface: 00026 C --------- 00027 C *CALL* *filling(pfild, pworka, pworkb, pworkc, klon, klat, 00028 C kmask, kmesh, kunit, cdfic, cdmet)* 00029 C 00030 C Input: 00031 C ----- 00032 C pfild : field to be completed (real 2D) 00033 C plon : grid longitude array (real 2D) 00034 C plat : grid latitude array (real 2D) 00035 C klon : number of longitude (integer) 00036 C klat : number of latitude (integer) 00037 C kmask : mask array (integer 2D) 00038 C kmesh : overlapped neighbors array (integer 2D) 00039 C kunit : climatological data filename logical unit (integer) 00040 C cdfic : climatological data filename (character) 00041 C cdmet : filling method (character) 00042 C 00043 C Output: 00044 C ------ 00045 C pfild : field completed (real 2D) 00046 C 00047 C Workspace: 00048 C --------- 00049 C These are work arrays passed as arguments: 00050 C pworka, pworkb, pworkc 00051 C 00052 C Externals: 00053 C --------- 00054 C None 00055 C 00056 C Reference: 00057 C --------- 00058 C See OASIS manual (1995) 00059 C 00060 C History: 00061 C ------- 00062 C Version Programmer Date Description 00063 C ------- ---------- ---- ----------- 00064 C 1.0 L. Terray 94/01/01 created 00065 C 1.1 L. Terray 94/10/01 modified: lecture of climagrd 00066 C 2.0beta L. Terray 95/12/10 modified: new structure 00067 C 2.0 L. Terray 96/02/01 modified: change on rewind() 00068 C 2.1 L. Terray 96/09/09 modified: flux correction term 00069 C 2.3 S. Valcke 99/04/30 added: printing levels 00070 C 00071 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00072 C 00073 C* ----------------Include files and USE of modules--------------------------- 00074 C 00075 USE mod_unit 00076 USE mod_parameter 00077 USE mod_analysis 00078 USE mod_coast 00079 USE mod_calendar 00080 USE mod_smooth 00081 USE mod_printing 00082 C 00083 C* ---------------------------- Argument declarations---------------------- 00084 C 00085 REAL (kind=ip_realwp_p) pfild(klon, klat), plon(klon, klat), 00086 $ plat(klon, klat) 00087 REAL (kind=ip_realwp_p) pworka(klon, klat), pworkb(klon, klat) 00088 REAL (kind=ip_realwp_p) pworkc(klon, klat) 00089 INTEGER (kind=ip_intwp_p) kmask(klon,klat), kmesh(klon, klat) 00090 CHARACTER*8 cdmet, cdfic 00091 C 00092 C* ---------------------------- Local declarations ---------------------- 00093 C 00094 CHARACTER*3 clsmooth, clfield 00095 CHARACTER*2 clclim 00096 C 00097 C* ---------------------------- Poema verses ---------------------------- 00098 C 00099 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00100 C 00101 C* 1. Initialization 00102 C -------------- 00103 C 00104 IF (nlogprt .GE. 2) THEN 00105 WRITE (UNIT = nulou,FMT = *) ' ' 00106 WRITE (UNIT = nulou,FMT = *) 00107 $ ' ROUTINE filling - Level 3' 00108 WRITE (UNIT = nulou,FMT = *) 00109 $ ' *************** *******' 00110 WRITE (UNIT = nulou,FMT = *) ' ' 00111 WRITE (UNIT = nulou,FMT = *) 00112 $ ' Filling and smoothing of field with climatology' 00113 WRITE (UNIT = nulou,FMT = *) ' ' 00114 CALL FLUSH(nulou) 00115 ENDIF 00116 C 00117 C* initialize error flag for error routine 00118 C 00119 iflag = 0 00120 C 00121 C* Get smoothing and time info 00122 C 00123 clsmooth = cdmet(1:3) 00124 clfield = cdmet(4:6) 00125 clclim = cdmet(7:8) 00126 C 00127 C* Get congelation temperature for SST and SIE treatment. 00128 C 00129 zmerg = -1.93 00130 C 00131 C 00132 C* 2. Getting climatological or interannual field 00133 C ------------------------------------------- 00134 C 00135 C* Monthly climatology 00136 C 00137 IF (nlogprt .GE. 2) THEN 00138 CALL prcout 00139 $ ('Complementary data file is file =', cdfic, 1) 00140 ENDIF 00141 IF (clclim .EQ. 'SE' .OR. clclim .EQ. 'MO') THEN 00142 IF (nlogprt .GE. 2) THEN 00143 WRITE (UNIT = nulou,FMT = *) ' ' 00144 ENDIF 00145 C 00146 C* Monthly climatology 00147 C 00148 IF (clclim .EQ. 'SE') THEN 00149 IF (nlogprt .GE. 2) THEN 00150 CALL prtout 00151 $ ( 00152 'Monthly records skipped in climatology file $ nsrec =', nsrec, 1) 00153 ENDIF 00154 C 00155 C* Read the right records 00156 C 00157 REWIND kunit 00158 IF (nsrec .eq. 1) THEN 00159 READ (kunit) pworka 00160 READ (kunit) pworkb 00161 ELSEIF (nsrec .gt. 1 .and. nsrec .lt. 12) THEN 00162 DO 210 ji = 1, nsrec - 1 00163 READ (kunit) pworka 00164 210 CONTINUE 00165 READ (kunit) pworka 00166 READ (kunit) pworkb 00167 ELSEIF (nsrec .eq. 12 .or. nsrec .eq. 0) THEN 00168 DO 220 ji = 1, 11 00169 READ (kunit) pworka 00170 220 CONTINUE 00171 READ (kunit) pworka 00172 REWIND kunit 00173 READ (kunit) pworkb 00174 ELSE 00175 WRITE (UNIT = nulou,FMT = *) ' ***WARNING***' 00176 WRITE (UNIT = nulou,FMT = *) 00177 $ ' ===>>> : wrong value for nsrec ' 00178 WRITE (UNIT = nulou,FMT = *) 00179 $ ' ====== ===== ' 00180 WRITE (UNIT = nulou,FMT = *) ' ' 00181 WRITE (UNIT = nulou,FMT = *) 00182 $ ' Variable nsrec = ',nsrec 00183 WRITE (UNIT = nulou,FMT = *) 00184 $ ' We STOP !!! Check value for variable nsrec ' 00185 CALL HALTE ('STOP in filling') 00186 ENDIF 00187 C 00188 C* Monthly interannual data 00189 C 00190 ELSE IF(clclim .EQ. 'MO') THEN 00191 IF (nlogprt .GE. 2) THEN 00192 CALL prtout 00193 $ ( 00194 'Monthly records skipped in interannual file $ nmrec =', nmrec, 1) 00195 ENDIF 00196 C 00197 C* Read the right records 00198 C 00199 REWIND kunit 00200 IF (nmrec .eq. 0) THEN 00201 READ (kunit) pworka 00202 READ (kunit) pworkb 00203 ELSEIF (nmrec .ge. 1) THEN 00204 DO 230 ji = 1, nmrec 00205 READ (kunit) pworka 00206 230 CONTINUE 00207 READ (kunit) pworka 00208 READ (kunit) pworkb 00209 ENDIF 00210 ENDIF 00211 C 00212 C* Get all coefficients for time linear interpolation 00213 C 00214 zpoid1 = FLOAT(njtwo - njnow) / FLOAT(njtwo - njone) 00215 zpoid2 = 1.0 - zpoid1 00216 zpoib1 = FLOAT(ndtwo - njnow) / FLOAT(ndtwo - ndone) 00217 zpoib2 = 1.0 - zpoib1 00218 C 00219 C* Get climatological field for current day 00220 C 00221 C* Specific treatment for SST and SIE due to freezing temperature 00222 C In both cases, one reads the temperature file and looks at 00223 C sea-ice transitions. One makes the assumption that the ice is 00224 C formed or disappears on the first of the month. 00225 C 00226 DO 240 jj = 1, klat 00227 DO 250 ji = 1, klon 00228 C 00229 C* General case 00230 C 00231 pworkc(ji,jj) = zpoid1 * pworka(ji,jj) + 00232 $ zpoid2 * pworkb(ji,jj) 00233 C 00234 C* If field is SST 00235 C 00236 IF (clfield .EQ. 'SST') THEN 00237 C 00238 C* Sea Ice to Sea 00239 C 00240 IF (pworka(ji,jj) .LE. zmerg .AND. 00241 $ pworkb(ji,jj) .GT. zmerg) THEN 00242 IF (njnow .GT. 15) THEN 00243 CONTINUE 00244 ELSE 00245 pworkc(ji,jj) = zmerg * zpoib1 + 00246 $ pworkb(ji,jj) * zpoib2 00247 ENDIF 00248 ENDIF 00249 C 00250 C* Sea to Sea Ice 00251 C 00252 IF (pworka(ji,jj) .GT. zmerg .AND. 00253 $ pworkb(ji,jj) .LE. zmerg) THEN 00254 IF (njnow .GT. 15) THEN 00255 pworkc(ji,jj) = pworka(ji,jj) * zpoib1 + 00256 $ zmerg * zpoib2 00257 ENDIF 00258 ENDIF 00259 C 00260 C* Sea to Sea 00261 C 00262 IF (pworka(ji,jj) .GT. zmerg .AND. 00263 $ pworkb(ji,jj) .GT. zmerg) THEN 00264 pworkc(ji,jj) = zpoid1 * pworka(ji,jj) + 00265 $ zpoid2 * pworkb(ji,jj) 00266 ENDIF 00267 ENDIF 00268 C 00269 C* If field is SIE 00270 C 00271 IF (clfield .EQ. 'SIE') THEN 00272 C 00273 C* Sea Ice to Sea Ice 00274 C 00275 IF (pworka(ji,jj) .LE. zmerg .AND. 00276 $ pworkb(ji,jj) .LE. zmerg) THEN 00277 pworkc(ji,jj) = 1.0 00278 ENDIF 00279 C 00280 C* Sea Ice to Sea 00281 C 00282 IF (pworka(ji,jj) .LE. zmerg .AND. 00283 $ pworkb(ji,jj) .GT. zmerg) THEN 00284 IF (njnow .gt. 15) THEN 00285 pworkc(ji,jj) = 1.0 00286 ELSE 00287 pworkc(ji,jj) = 0.0 00288 ENDIF 00289 ENDIF 00290 C 00291 C* Sea to Sea Ice 00292 C 00293 IF (pworka(ji,jj) .GT. zmerg .AND. 00294 $ pworkb(ji,jj) .LE. zmerg) THEN 00295 IF (njnow .gt. 15) THEN 00296 pworkc(ji,jj) = 0.0 00297 ELSEIF (njnow .eq. 1) THEN 00298 pworkc(ji,jj) = 2.0 00299 ELSE 00300 pworkc(ji,jj) = 1.0 00301 ENDIF 00302 ENDIF 00303 C 00304 C* Sea to Sea 00305 C 00306 IF (pworka(ji,jj) .GT. zmerg .AND. 00307 $ pworkb(ji,jj) .GT. zmerg) THEN 00308 pworkc(ji,jj) = 0.0 00309 ENDIF 00310 ENDIF 00311 250 CONTINUE 00312 240 CONTINUE 00313 C 00314 C* Annual climatology 00315 C 00316 ELSE IF (clclim .EQ. 'AN') THEN 00317 REWIND kunit 00318 READ (kunit) pworkc 00319 ELSE 00320 CALL prcout ( 00321 'This type of climatology cannot be used - $ clim = ', clclim, 1) 00322 ENDIF 00323 C 00324 C* Blending and smoothing 00325 C 00326 C* Storing of climatology field interpolated in intermediate array pworka 00327 C 00328 DO 260 jj = 1, klat 00329 DO 270 ji = 1, klon 00330 pworka(ji,jj) = pworkc(ji,jj) 00331 270 CONTINUE 00332 260 CONTINUE 00333 C* SST case 00334 IF (clfield .EQ. 'SST') THEN 00335 C 00336 C* Put interpolated values in intermediate array if: 00337 C - the point is an atmosphere sea point 00338 C - there is at least one underlying ocean sea point 00339 C 00340 DO 280 jj = 1, klat 00341 DO 290 ji = 1, klon 00342 IF (kmesh(ji,jj) .NE. 0 .AND. kmask(ji,jj) .EQ. 0) THEN 00343 pworka(ji,jj) = pfild(ji,jj) 00344 ENDIF 00345 290 CONTINUE 00346 280 CONTINUE 00347 C 00348 C* If needed, do coast mismatch correction 00349 C --> Special treatment for coastal points where there is a 00350 C --> a mismatch between ocean and atmosphere land sea masks 00351 C --> The array npcoast is calculated in routine coasts.f 00352 C --> icoor and jcoor are externally defined functions 00353 C 00354 IF(nfcoast .EQ. 1) THEN 00355 C* If first iteration, initialize coast mismatch data 00356 IF (lcoast) THEN 00357 itmsq = 1 00358 CALL coasts (plon, plat, kmask, itmsq, 00359 $ klon, klat, kmesh) 00360 C* Switch off flag 00361 lcoast = .FALSE. 00362 ENDIF 00363 IF (nlogprt .GE. 2) THEN 00364 WRITE (UNIT = nulan,FMT = *) ' ' 00365 WRITE (UNIT = nulan,FMT = *) 00366 $ ' * Model run date * ' 00367 WRITE (UNIT = nulan,FMT = *) 00368 $ ' * -------------- * ' 00369 WRITE (UNIT = nulan,FMT = *) ' ' 00370 WRITE (UNIT = nulan,FMT = *) 00371 $ ' Year --->>> ',nanow 00372 WRITE (UNIT = nulan,FMT = *) 00373 $ ' Month --->>> ',nmnow 00374 WRITE (UNIT = nulan,FMT = *) 00375 $ ' Day --->>> ',njnow 00376 WRITE (UNIT = nulan,FMT = *) ' ' 00377 WRITE (UNIT = nulan,FMT = *) 00378 $ ' Correct for coast mismatch ' 00379 WRITE (UNIT = nulan,FMT = *) 00380 $ ' ************************** ' 00381 WRITE (UNIT = nulan,FMT = *) ' ' 00382 ENDIF 00383 DO 291 jc = 1, ncoast 00384 iic = icoor(npcoast(jc,1), klon) 00385 ijc = jcoor(npcoast(jc,1), klon) 00386 ztmptot = 0. 00387 DO 292 jn = 3, 2 + npcoast(jc,2) 00388 iiic = icoor(npcoast(jc,jn), klon) 00389 ijjc = jcoor(npcoast(jc,jn), klon) 00390 ztmptot = ztmptot + pfild(iiic,ijjc) 00391 292 CONTINUE 00392 C 00393 C* Print old and new values 00394 C 00395 IF (nlogprt .GE. 2) THEN 00396 WRITE (UNIT = nulan,FMT = *) 00397 $ ' Resetting SST at point (i,j)= ',iic,ijc 00398 WRITE (UNIT = nulan,FMT = *) 00399 $ ' Old SST value : SST = ', pfild(iic,ijc) 00400 WRITE (UNIT = nulan,FMT = *) 00401 $ ' Clim SST value : SST = ', pworka(iic,ijc) 00402 ENDIF 00403 C 00404 C* Get new value 00405 C 00406 pfild(iic,ijc) = ztmptot / float(npcoast(jc,2)) 00407 IF (nlogprt .GE. 2) THEN 00408 WRITE (UNIT = nulan,FMT = *) 00409 $ ' New SST value : SST = ', pfild(iic,ijc) 00410 ENDIF 00411 C 00412 C* Copy new value into intermediate array 00413 C 00414 pworka(iic,ijc) = pfild(iic,ijc) 00415 291 CONTINUE 00416 ENDIF 00417 C 00418 C* SIE case : we assign the new field to climatology i.e we do nothing 00419 C 00420 ELSE IF (clfield .EQ. 'SIE') THEN 00421 CONTINUE 00422 C 00423 C* General case: put interpolated value if there is an 00424 C underlying grid point (whether land or sea) 00425 C 00426 ELSE 00427 DO 293 jj = 1, klat 00428 DO 294 ji = 1, klon 00429 IF (kmesh(ji,jj) .NE. 0) THEN 00430 pworka(ji,jj) = pfild(ji,jj) 00431 ENDIF 00432 294 CONTINUE 00433 293 CONTINUE 00434 ENDIF 00435 C 00436 C 00437 C* 3. Do the smoothing on borders (only for SST) 00438 C ---------------------------- 00439 C 00440 IF (clfield .EQ. 'SST' .AND. clsmooth .EQ. 'SMO') THEN 00441 C 00442 C* First South and North borders 00443 C 00444 DO 310 jk = nsltb, nslte 00445 DO 320 ji = 1, klon 00446 IF (kmesh(ji,jk) .NE. 0 .AND. kmask(ji,jk) .EQ. 0) THEN 00447 pworka(ji,jk) = qalfa * (jk - nsltb) * pfild(ji,jk) 00448 $ + (1.-qalfa * (jk - nsltb)) * pworkc(ji,jk) 00449 ENDIF 00450 320 CONTINUE 00451 310 CONTINUE 00452 DO 330 jk = nnltb, nnlte, -1 00453 DO 340 ji = 1, klon 00454 IF (kmesh(ji,jk) .NE. 0 .AND. kmask(ji,jk) .EQ. 0) THEN 00455 pworka(ji,jk) = qalfa * (nnltb - jk) * pfild(ji,jk) 00456 $ + (1.-qalfa * (nnltb - jk)) * pworkc(ji,jk) 00457 ENDIF 00458 340 CONTINUE 00459 330 CONTINUE 00460 C 00461 C* Western border 00462 C 00463 DO 350 jj = nslte + 1, nnlte - 1 00464 DO 360 ji = 1, nwlgmx 00465 IF (kmesh(ji,jj) .EQ. 0 .AND. 00466 $ kmesh(ji + 1,jj) .NE. 0) THEN 00467 DO 370 jk = 1, nliss 00468 pworka(ji + jk - 1,jj) = qbeta*(jk - 1) * 00469 $ pfild(ji + jk - 1,jj) + 00470 $ pworkc(ji + jk - 1,jj) * 00471 $ (1.-qbeta*(jk - 1)) 00472 370 CONTINUE 00473 ENDIF 00474 360 CONTINUE 00475 350 CONTINUE 00476 C 00477 C* Eastern border 00478 C 00479 DO 380 jj = nslte + 1, nnlte - 1 00480 DO 390 ji = klon, nelgmx, -1 00481 IF (kmesh(ji,jj) .EQ. 0 .AND. 00482 $ kmesh(ji - 1,jj) .NE. 0) THEN 00483 DO 391 jk = 1, nliss 00484 pworka(ji - jk + 1,jj) = qbeta*(jk - 1) * 00485 $ pfild(ji - jk + 1,jj) + 00486 $ pworkc(ji - jk + 1,jj) * 00487 $ (1.-qbeta*(jk - 1)) 00488 391 CONTINUE 00489 ENDIF 00490 390 CONTINUE 00491 380 CONTINUE 00492 C 00493 C* Calculate and write flux correction term cfldcor on unit nlucor 00494 C 00495 CALL prcout('Calculate flux correction term', cfldcor, 1) 00496 CALL szero(pworkb,klon*klat) 00497 DO 392 jj = 1, klat 00498 DO 393 ji = 1, klon 00499 C 00500 C* The correction term is estimated only for atmosphere sea points with 00501 C underlying ocean sea points 00502 C 00503 IF (kmesh(ji,jj) .NE. 0 .AND. kmask(ji,jj) .EQ. 0) THEN 00504 pworkb(ji,jj) = pfild(ji,jj) - pworka(ji,jj) 00505 ENDIF 00506 393 CONTINUE 00507 392 CONTINUE 00508 C 00509 C* Rewind unit nlucor and write 00510 C 00511 IF (nlogprt .GE. 2) THEN 00512 CALL prtout('Write it on logical unit ', nlucor, 2) 00513 ENDIF 00514 REWIND nlucor 00515 CALL locwrite(cfldcor, pworkb, klon*klat, nlucor, iflag) 00516 IF (iflag .NE. 0) THEN 00517 CALL prcout 00518 $ ('WARNING: problem in writing field', 00519 $ cfldcor, 1) 00520 CALL prtout 00521 $ ('Error in writing logical unit', nlucor, 2) 00522 CALL HALTE('STOP in filling') 00523 ENDIF 00524 ENDIF 00525 C 00526 C* Assign intermediate array pworka values to final array pfild 00527 C 00528 DO 394 jj = 1, klat 00529 DO 395 ji = 1, klon 00530 pfild(ji,jj) = pworka(ji,jj) 00531 395 CONTINUE 00532 394 CONTINUE 00533 C 00534 C 00535 C* 4. End of routine 00536 C -------------- 00537 C 00538 IF (nlogprt .GE. 2) THEN 00539 WRITE (UNIT = nulou,FMT = *) ' ' 00540 WRITE (UNIT = nulou,FMT = *) 00541 $ ' --------- End of ROUTINE filling ---------' 00542 CALL FLUSH (nulou) 00543 ENDIF 00544 RETURN 00545 END
1.7.4 
