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