|
Oasis3 4.0.2
|
|
Oasis3 4.0.2
|
00001 SUBROUTINE namset (px1, py1, kmsk1, kvmsk1, kngx1, kngy1, 00002 $ cdper1, kper1, 00003 $ px2, py2, kmsk2, kvmsk2, kngx2, kngy2, 00004 $ cdper2, kper2, 00005 $ cdtyp, 00006 $ pr1to2, k1to2, kw1to2, psgr1, psgr2, 00007 $ krdwt, kwlun, kvmsz2, 00008 $ kmskz2, knumber) 00009 C**** 00010 C ***************************** 00011 C * OASIS ROUTINE - LEVEL 3 * 00012 C * ------------- ------- * 00013 C ***************************** 00014 C 00015 C**** *namset* - Initializations for the mesh naive method 00016 C 00017 C Purpose: 00018 C ------- 00019 C Read or calculate the interpolation weights used in Anaism 00020 C 00021 C** Interface: 00022 C --------- 00023 C *CALL* *namset(px1, py1, kmsk1, kvmsk1, kngx1, kngy1, cdper1, kper1, 00024 C px2, py2, kmsk2, kvmsk2, kngx2, kngy2, cdper2, kper2, 00025 C pr1to2, k1to2, kw1to2, psgr1, psgr2, 00026 C krdwt, kwlun, kvmsz2, 00027 C kmskz2, knumber)* 00028 C Input: 00029 C ----- 00030 C kngx1 : number of longitudes for source grid 00031 C kngy1 : number of latitudes for source grid 00032 C px1 : longitudes for source grid (real 2D) 00033 C py1 : latitudes for source grid (real 2D) 00034 C kmsk1 : the mask for source grid (integer 2D) 00035 C kvmsk1 : the value of the mask for source grid 00036 C cdper1 : source grid periodicity 00037 C kper1 : number of overlapped points for source grid 00038 C kngx2 : number of longitudes for target grid 00039 C kngy2 : number of latitudes for target grid 00040 C px2 : longitudes for target grid (real 2D) 00041 C py2 : latitudes for target grid (real 2D) 00042 C kmsk2 : the mask of target grid (integer 2D) 00043 C kvmsk2 : the value of the mask for target grid 00044 C cdper2 : target grid periodicity 00045 C kper2 : number of overlapped points for target grid 00046 C cdtyp : type of source grid 00047 C kw1to2 : maximum number of overlapped neighbors 00048 C krdwt : read/write flag for the weights 00049 C kwlun : logical unit for the weights 00050 C kvmsz2 : mask value for array kmskz2 00051 C knumber : flag to identify appropriate Anaism dataset 00052 C 00053 C Output: 00054 C ------ 00055 C pr1to2 : weights for Anaism interpolation (real 2D) 00056 C k1to2 : source grid neighbors adresses (integer 2D) 00057 C kmskz2 : number of source grid neighbors (integer 2D) 00058 C psgr1 : source grid square surfaces (real 2D) 00059 C psgr2 : target grid square surfaces (real 2D) 00060 C 00061 C 00062 C Workspace: 00063 C --------- 00064 C None 00065 C 00066 C External: 00067 C -------- 00068 C pcssph, ssumr, grstat, pmrhal, locwrint, locwrite, locread, locrint, 00069 C icoor, jcoor 00070 C 00071 C References: 00072 C ---------- 00073 C O. Thual, Simple ocean-atmosphere interpolation. 00074 C Part A: The method, EPICOA 0629 (1992) 00075 C Part B: Software implementation, EPICOA 0630 (1992) 00076 C See also OASIS manual (1995) 00077 C 00078 C History: 00079 C ------- 00080 C Version Programmer Date Description 00081 C ------- ---------- ---- ----------- 00082 C 1.1 O. Thual 93/04/15 created 00083 C 2.0 L. Terray 95/10/01 modified: new structure 00084 C 2.3 S. Valcke 99/04/30 added: printing levels 00085 C 2.3 L. Terray 99/09/15 changed periodicity variables 00086 C 00087 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00088 C 00089 C* ---------------------------- Include files --------------------------- 00090 C 00091 USE mod_unit 00092 USE mod_printing 00093 C 00094 C* ---------------------------- Argument declarations ------------------- 00095 C 00096 REAL (kind=ip_realwp_p) px1(kngx1,kngy1), py1(kngx1,kngy1), 00097 $ psgr1(kngx1,kngy1) 00098 REAL (kind=ip_realwp_p) px2(kngx2,kngy2), py2(kngx2,kngy2), 00099 $ psgr2(kngx2,kngy2) 00100 REAL (kind=ip_realwp_p) pr1to2(kw1to2,kngx2*kngy2) 00101 INTEGER (kind=ip_intwp_p) kmsk1(kngx1,kngy1), kmsk2(kngx2,kngy2) 00102 INTEGER (kind=ip_intwp_p) k1to2(kw1to2,kngx2*kngy2), 00103 $ kmskz2(kngx2,kngy2) 00104 CHARACTER*8 cdper1, cdper2 00105 CHARACTER*1 cdtyp 00106 C 00107 C* ---------------------------- Local declarations ---------------------- 00108 C 00109 CHARACTER*8 cladress, clweight, cloverlp 00110 CHARACTER*1 cltyp 00111 C 00112 C* ---------------------------- Poema verses ---------------------------- 00113 C 00114 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00115 C 00116 C* 1. Initializations 00117 C --------------- 00118 C 00119 IF (nlogprt .GE. 2) THEN 00120 WRITE (UNIT = nulou,FMT = *) ' ' 00121 WRITE (UNIT = nulou,FMT = *) ' ' 00122 WRITE (UNIT = nulou,FMT = *) 00123 $ ' ROUTINE namset - Level 3' 00124 WRITE (UNIT = nulou,FMT = *) 00125 $ ' ************** *******' 00126 WRITE (UNIT = nulou,FMT = *) ' ' 00127 WRITE (UNIT = nulou,FMT = *) 00128 $ ' Set up ANAIS-MESH interpolation' 00129 WRITE (UNIT = nulou,FMT = *) ' ' 00130 WRITE (UNIT = nulou,FMT = *) ' ' 00131 CALL FLUSH(nulou) 00132 ENDIF 00133 C 00134 C* Define global dimensions + other local variables 00135 C 00136 ing1 = kngx1 * kngy1 00137 ing2 = kngx2 * kngy2 00138 ivmsk1 = kvmsk1 00139 ivmsk2 = kvmsk2 00140 iflag = 0 00141 C 00142 C* Define character strings to locate needed data 00143 C 00144 WRITE(clweight,'(''WEIGHTS'',I1)') knumber 00145 WRITE(cladress,'(''ADRESSE'',I1)') knumber 00146 WRITE(cloverlp,'(''OVERLAP'',I1)') knumber 00147 C 00148 C 00149 C* 2. Statistics of source grid 00150 C ------------------------- 00151 C 00152 cltyp = cdtyp 00153 C* The following routines calculate some interesting info about the grids 00154 C 00155 C* Calculate surface elements (assume spherical and periodic grid) 00156 C 00157 CALL pcssph (px1, py1, psgr1, kngx1, kngy1) 00158 zsum = ssumr (psgr1, ing1) 00159 C 00160 C* Printing on ANAIS output file 00161 C 00162 IF (nlogprt .GE. 2) THEN 00163 WRITE (UNIT = nulan,FMT = *) ' ANAIS output FILE ' 00164 WRITE (UNIT = nulan,FMT = *) ' ***** ****** **** ' 00165 WRITE (UNIT = nulan,FMT = *) ' ' 00166 WRITE (UNIT = nulan,FMT = *) ' Routine namset ' 00167 WRITE (UNIT = nulan,FMT = *) ' -------------- ' 00168 WRITE (UNIT = nulan,FMT = *) ' ' 00169 WRITE (UNIT = nulan,FMT = *) 00170 $ ' Some statistics about the source grid ' 00171 WRITE (UNIT = nulan,FMT = *) 00172 $ ' ************************************* ' 00173 WRITE (UNIT = nulan,FMT = *) ' ' 00174 WRITE (UNIT = nulan,FMT = *) 00175 $ ' Sum of grid square surfaces = ', zsum 00176 WRITE (UNIT = nulan,FMT = *) 00177 $ ' It must be equal to 4 x PI = ', 16. * atan(1.) 00178 WRITE (UNIT = nulan,FMT = *) ' ' 00179 ENDIF 00180 C 00181 C* Calculate for ocean points (unmasked points) the following data : 00182 C ************ 00183 C - total surface 00184 C - average inter-point distance 00185 C - inverse of the sum of surface elements squared 00186 C 00187 CALL grstat (px1, py1, psgr1, kmsk1, kngx1, kngy1, cltyp, 00188 $ zamsh1, zstot1, zhhi1, ivmsk1) 00189 C 00190 C* Printing 00191 C 00192 IF (nlogprt .GE. 2) THEN 00193 WRITE (UNIT = nulan,FMT = *) 00194 $ ' Some info for ocean points only ' 00195 WRITE (UNIT = nulan,FMT = *) 00196 $ ' ******************************* ' 00197 WRITE (UNIT = nulan,FMT = *) ' ' 00198 WRITE (UNIT = nulan,FMT = *) 00199 $ ' Average mesh distance = ', sqrt(zamsh1) 00200 WRITE (UNIT = nulan,FMT = *) ' Total surface = ', zstot1 00201 WRITE (UNIT = nulan,FMT = *) 00202 $ ' Inverse of sum surf. elts. squared = ', zhhi1 00203 ENDIF 00204 C 00205 C 00206 C* 3. Statistics of target grid 00207 C ------------------------- 00208 C 00209 C* The following routines calculate some interesting info about the grids 00210 C 00211 C* Calculate surface elements (assume spherical and periodic grid) 00212 C 00213 CALL pcssph (px2, py2, psgr2, kngx2, kngy2) 00214 zsum = ssumr (psgr2, ing2) 00215 C 00216 C* Printing 00217 C 00218 IF (nlogprt .GE. 2) THEN 00219 WRITE (UNIT = nulan,FMT = *) ' ' 00220 WRITE (UNIT = nulan,FMT = *) 00221 $ ' Some statistics about the target grid ' 00222 WRITE (UNIT = nulan,FMT = *) 00223 $ ' ************************************* ' 00224 WRITE (UNIT = nulan,FMT = *) ' ' 00225 WRITE (UNIT = nulan,FMT = *) 00226 $ ' Sum of grid square surfaces = ', zsum 00227 WRITE (UNIT = nulan,FMT = *) 00228 $ ' It must be equal to 4 x PI = ', 00229 $ 16. * atan(1.) 00230 WRITE (UNIT = nulan,FMT = *) ' ' 00231 ENDIF 00232 C 00233 C* Calculate for ocean points (unmasked points) the following data : 00234 C ************ 00235 C - total surface 00236 C - average inter-point distance 00237 C - inverse of the sum of surface elements squared 00238 C 00239 CALL grstat (px2, py2, psgr2, kmsk2, kngx2, kngy2, cltyp, 00240 $ zamsh2, zstot2, zhhi2, ivmsk2) 00241 C 00242 C* Printing 00243 C 00244 IF (nlogprt .GE. 2) THEN 00245 WRITE (UNIT = nulan,FMT = *) 00246 $ ' Some info for ocean points only ' 00247 WRITE (UNIT = nulan,FMT = *) 00248 $ ' ******************************* ' 00249 WRITE (UNIT = nulan,FMT = *) ' ' 00250 WRITE (UNIT = nulan,FMT = *) 00251 $ ' Average mesh distance = ', sqrt(zamsh2) 00252 WRITE (UNIT = nulan,FMT = *) ' Total surface = ', zstot2 00253 WRITE (UNIT = nulan,FMT = *) 00254 $ ' Inverse of sum surf. elts. squared = ', zhhi2 00255 ENDIF 00256 C 00257 C 00258 C* 4. Preparation of the neighbors search 00259 C ----------------------------------- 00260 C 00261 C* Consistency checking 00262 C 00263 IF (kw1to2 .GT. ing2) THEN 00264 WRITE (UNIT = nulou,FMT = *) ' WARNING: kw1to2 .gt. ing2 ' 00265 WRITE (UNIT = nulou,FMT = *) ' ******* ------ ---- ' 00266 WRITE (UNIT = nulou,FMT = *) 00267 $ ' kw1to2 = ',kw1to2,' ing2 = ',ing2 00268 ENDIF 00269 C 00270 C 00271 C* 5. Weights determination 00272 C --------------------- 00273 C 00274 C* Initialization for read write options 00275 C 00276 IF (krdwt .EQ. 1) THEN 00277 C 00278 C* Writing 00279 C ------- 00280 C 00281 C* Calculate and write weights after checking numbers 00282 C 00283 C* Calculates weights 00284 C 00285 CALL pmrhal (pr1to2, k1to2, kw1to2, 00286 $ px1, py1, kmsk1, kngx1, kngy1, cdper1, kper1, 00287 $ px2, py2, kmsk2, kngx2, kngy2, cdper2, kper2, 00288 $ ivmsk1, ivmsk2, kmskz2, kvmsz2) 00289 C 00290 C* Write weights + other useful stuff 00291 C 00292 C* - Write main locator string + weights 00293 C 00294 CALL locwrite (clweight, pr1to2, kw1to2*kngx2*kngy2, 00295 $ kwlun, iflag) 00296 IF (iflag .NE. 0) THEN 00297 WRITE (UNIT = nulou,FMT = *) 00298 $ 'Problem in writing on UNIT = ', kwlun 00299 WRITE (UNIT = nulou,FMT = *) 00300 $ 'String locator is = ', clweight 00301 CALL HALTE ('Stop in namset') 00302 ENDIF 00303 C 00304 C* - Write overlapped neighbors adresses 00305 C 00306 CALL locwrint (cladress, k1to2, kw1to2*kngx2*kngy2, 00307 $ kwlun, iflag) 00308 IF (iflag .NE. 0) THEN 00309 WRITE (UNIT = nulou,FMT = *) 00310 $ 'Problem in writing on UNIT = ', kwlun 00311 WRITE (UNIT = nulou,FMT = *) 00312 $ 'String locator is = ', cladress 00313 CALL HALTE ('Stop in namset') 00314 ENDIF 00315 C 00316 C* - Write overlapped neighbors numbers 00317 C 00318 CALL locwrint (cloverlp, kmskz2, kngx2*kngy2, 00319 $ kwlun, iflag) 00320 IF (iflag .NE. 0) THEN 00321 WRITE (UNIT = nulou,FMT = *) 00322 $ 'Problem in writing on UNIT = ', kwlun 00323 WRITE (UNIT = nulou,FMT = *) 00324 $ 'String locator is = ', cloverlp 00325 CALL HALTE ('Stop in namset') 00326 ENDIF 00327 C 00328 C* - Write checkings numbers (grid parameters) 00329 C 00330 WRITE(kwlun) kngx1, kngy1, kngx2, kngy2 00331 C 00332 C* Job is done 00333 C 00334 IF (nlogprt .GE. 2) THEN 00335 CALL prtout('Wrote weights on unit = ', kwlun, 1) 00336 ENDIF 00337 C 00338 C* Reading 00339 C ------- 00340 C 00341 ELSE 00342 C 00343 C 00344 C* Reads weights and checks 00345 C 00346 C* - Read main locator string + weights 00347 CALL locread (clweight, pr1to2, kw1to2*kngx2*kngy2, 00348 $ kwlun, iflag) 00349 IF (iflag .NE. 0) THEN 00350 WRITE (UNIT = nulou,FMT = *) 00351 $ 'Problem in reading on UNIT = ', kwlun 00352 WRITE (UNIT = nulou,FMT = *) 00353 $ 'String locator is = ', clweight 00354 CALL HALTE ('Stop in namset') 00355 ENDIF 00356 C 00357 C* - Read overlapped neighbors adresses 00358 C 00359 CALL locrint (cladress, k1to2, kw1to2*kngx2*kngy2, 00360 $ kwlun, iflag) 00361 IF (iflag .NE. 0) THEN 00362 WRITE (UNIT = nulou,FMT = *) 00363 $ 'Problem in reading on UNIT = ', kwlun 00364 WRITE (UNIT = nulou,FMT = *) 00365 $ 'String locator is = ', cladress 00366 CALL HALTE ('Stop in namset') 00367 ENDIF 00368 C 00369 C* - Read overlapped neighbors numbers 00370 C 00371 CALL locrint (cloverlp, kmskz2, kngx2*kngy2, 00372 $ kwlun, iflag) 00373 IF (iflag .NE. 0) THEN 00374 WRITE (UNIT = nulou,FMT = *) 00375 $ 'Problem in reading on UNIT = ', kwlun 00376 WRITE (UNIT = nulou,FMT = *) 00377 $ 'String locator is = ', cloverlp 00378 CALL HALTE ('Stop in namset') 00379 ENDIF 00380 C 00381 C* - Read checkings numbers (grid parameters) 00382 C 00383 READ(kwlun) ingx1, ingy1, ingx2, ingy2 00384 C 00385 C* Checks 00386 C 00387 IF (ingx1 .NE. kngx1 .OR. ingy1 .NE. kngy1 .OR. 00388 $ ingx2 .NE. kngx2 .OR. ingy2 .NE. kngy2) THEN 00389 WRITE (UNIT = nulou,FMT = *) 00390 $ ' Inconsistency in mweights file' 00391 WRITE (UNIT = nulou,FMT = *) 00392 $ 'ingx1 = ',ingx1,'kngx1 = ',kngx1 00393 WRITE (UNIT = nulou,FMT = *) 00394 $ 'ingy1 = ',ingy1,'kngy1 = ',kngy1 00395 WRITE (UNIT = nulou,FMT = *) 00396 $ 'ingx2 = ',ingx2,'kngx2 = ',kngx2 00397 WRITE (UNIT = nulou,FMT = *) 00398 $ 'ingy2 = ',ingy2,'kngy2 = ',kngy2 00399 CALL HALTE ('Stop in namset') 00400 ENDIF 00401 C 00402 C* Reading of weights done 00403 C 00404 IF (nlogprt .GE. 2) THEN 00405 CALL prtout 00406 $ ('Reading of weights done on unit = ', kwlun, 1) 00407 ENDIF 00408 ENDIF 00409 C 00410 C 00411 C* 6. Printing weights and adresses on ANAIS output file 00412 C -------------------------------------------------- 00413 C 00414 IF (nlogprt .GE. 2) THEN 00415 WRITE (UNIT = nulan,FMT = *) ' ' 00416 WRITE (UNIT = nulan,FMT = *) 00417 $ ' Print weights and adresses of overlapped neighbors ' 00418 WRITE (UNIT = nulan,FMT = *) 00419 $ ' ************************************************** ' 00420 WRITE (UNIT = nulan,FMT = *) ' ' 00421 DO 610 jj = 1, kngy2 00422 DO 620 ji = 1, kngx2 00423 WRITE (UNIT = nulan,FMT = 6010) ji, jj 00424 WRITE (UNIT = nulan,FMT = 6020) px2(ji,jj), py2(ji,jj) 00425 WRITE (UNIT = nulan,FMT = 6030) kmskz2(ji,jj) 00426 IF (kmskz2(ji,jj) .NE. 0) THEN 00427 IF (kmsk2(ji,jj) .EQ. 0) THEN 00428 DO 630 jn = 1, kmskz2(ji,jj) 00429 iind = ji + kngx2 * (jj-1) 00430 iadr = k1to2(jn,iind) 00431 ix = icoor(iadr,kngx1) 00432 iy = jcoor(iadr,kngx1) 00433 WRITE(UNIT = nulan,FMT = 6040) 00434 $ jn, pr1to2(jn,iind) 00435 WRITE(UNIT = nulan,FMT = 6050) 00436 $ ix, iy, px1(ix,iy), py1(ix,iy) 00437 630 CONTINUE 00438 ELSE 00439 WRITE(UNIT = nulan,FMT = 6060) 00440 ENDIF 00441 ENDIF 00442 620 CONTINUE 00443 610 CONTINUE 00444 CALL FLUSH(nulan) 00445 ENDIF 00446 C 00447 C* Formats 00448 C 00449 6010 FORMAT(/,' Target grid point -- i = ',I3,' j = ',I3) 00450 6020 FORMAT(' Point coordinates -- long = ',F9.4,' lat = ',F9.4) 00451 6030 FORMAT(3X,'Number of overlapped source grid squares = ',I3) 00452 6040 FORMAT(5X,' neighbor number = ',I3,' weight is = ',F6.4) 00453 6050 FORMAT(5X,' i = ',I3,' j = ',I3,' lon = ',F9.4,' lat = ',F9.4) 00454 6060 FORMAT(5X,' This is a continental point on the target grid ') 00455 C 00456 C 00457 C* 7. End of routine 00458 C -------------- 00459 C 00460 IF (nlogprt .GE. 2) THEN 00461 WRITE(UNIT = nulou,FMT = *) ' ' 00462 WRITE(UNIT = nulou,FMT = *) 00463 $ ' --------- End of ROUTINE namset ---------' 00464 CALL FLUSH(nulou) 00465 ENDIF 00466 RETURN 00467 END 00468 00469 00470 00471
1.7.4 
