|
Oasis3 4.0.2
|
|
Oasis3 4.0.2
|
00001 SUBROUTINE scriprmp_vector( & 00002 dst_arrayI_dstA, dst_arrayJ_dstB, & 00003 src_arrayI_srcA, src_arrayJ_srcB, & 00004 grd_name_srcA, grd_name_srcB, & 00005 grd_name_dstA, grd_name_dstB, & 00006 src_size, dst_size, & 00007 msk_srcA, msk_srcB, & 00008 msk_dstA, msk_dstB, & 00009 lon_srcA, lat_srcA, & 00010 lon_srcB, lat_srcB, & 00011 lon_dstA, lat_dstA, & 00012 lon_dstB, lat_dstB, & 00013 nlon_src, nlat_src, & 00014 nlon_dst, nlat_dst, & 00015 map_method, cdgrdtyp, & 00016 id_sper, id_tper, & 00017 cd_sper, cd_tper, & 00018 normalize_opt, order, & 00019 rst_type, n_srch_bins, & 00020 lextrapdone, lprojcart, & 00021 rl_varmul, id_scripvoi ) 00022 00023 !C Input: 00024 !C ----- 00025 !C src_arrayI_srcA : field on source grid defined as VECTOR_I 00026 !C src_arrayJ_srcB : field on source grid defined as VECTOR_J 00027 !C grd_name_srcA : name for the source grid where VECTOR_I is defined 00028 !C grd_name_srcB : name for the source grid where VECTOR_J is defined 00029 !C grd_name_dstA : name for the target grid where VECTOR_I is defined 00030 !C grd_name_dstB : name for the target grid where VECTOR_J is defined 00031 !C src_size : source grid size 00032 !C dst_size : target grid size 00033 !C msk_srcA : grid mask for srcA 00034 !C msk_srcB : grid mask for srcB 00035 !C msk_dstA : grid mask for dstA 00036 !C msk_dstB : grid mask for dstB 00037 !C lon_srcA : longitudes for grid srcA 00038 !C lon_srcB : longitudes for grid srcB 00039 !C lon_dstA : longitudes for grid dstA 00040 !C lon_dstB : longitudes for grid dstB 00041 !C lat_srcA : latitudes for grid srcA 00042 !C lat_srcB : latitudes for grid srcB 00043 !C lat_dstA : latitudes for grid dstA 00044 !C lat_dstB : latitudes for grid dstB 00045 !C nlon_src : number of longitudes for source grid 00046 !C nlat_src : number of latitudes for source grid 00047 !C nlon_dst : number of longitudes for target grid 00048 !C nlat_dst : number of latitudes for target grid 00049 !C map_method : remapping method 00050 !C cdgrdtyp : grid type for source grid 00051 !C id_sper : number of overlapping for source grid 00052 !C id_tper : number of overlapping for target grid 00053 !C cd_sper : source grid periodicity type 00054 !C cd_tper : target grid periodicity type 00055 !C normalize_opt : option for normalization 00056 !C order : order of conservative remapping 00057 !C rst_type : type of scrip search restriction 00058 !C n_srch_bins : number of seach bins 00059 !C lextrapdone : logical, true if EXTRAP done on field 00060 !C lprojcart : logical, true for rotation to cartesian refernetial 00061 !C rl_varmul : Gaussian variance (for GAUSWGT) 00062 !C id_scripvoi : number of neighbour for DISTWGT and GAUSWGT 00063 !C 00064 !C Output: 00065 !C ------ 00066 !C dst_arrayI_dstA : VECTOR_I field on target grid as defined in namcouple 00067 !C dst_arrayJ_dstB : VECTOR_J field on target grid as defined in namcouple 00068 !C 00069 !C Externals: 00070 !C --------- 00071 !C corners, scrip, gradient, gradient_bicubic, 00072 !C from module rotations : loc2spher, sher2car, car2spher, spher2loc 00073 !C from module vector : remap_vector_comps, write_src_array_spheric, 00074 !C write_dst_array_spheric, write_cartesian_components, 00075 !C check_points_at_poles 00076 !C 00077 !C History: 00078 !C -------- 00079 !C E. Rapaport 2004/03 Created 00080 !C J. Ghattas 2006/02 Rewritten 00081 ! 00082 ! 00083 ! Description : 00084 ! ------------- 00085 ! This routine interpolates a vector field, described as 2 scalar components I and J. 00086 ! The 2 components are described on different source grids and the resulting 2 00087 ! components can be calculated on 2 differents target grids. The source grids 00088 ! are named grd_srcA which is the original grid for the component I, and grd_srcB 00089 ! which is the original grid for component J. The target grid defined for the component 00090 ! I are dstA and for for the component J dstB. The two source grids have to be identical 00091 ! in number of cells and their mask must be the same; for different target grids the same 00092 ! criteria are imposed. 00093 ! 00094 ! Steps: 00095 ! ------ 00096 ! -> If the source grid differs a first interpolation will be done from grd_srcB 00097 ! towards grd_srcA scalarly for the 2 components 00098 ! -> Rotation towards spherical referentials on source grid srcA 00099 ! -> If desired rotation towards cartesian referentials 00100 ! -> Interpolation of the components scalarly from srcA towards target grid dstA 00101 ! -> Rotation towards local referentials. 00102 ! -> If target grids differ, the same procedure will be done from grd_srcA towards grd_dstB. 00103 ! 00104 !C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00105 !C* ---------------------------- Modules used ---------------------------- 00106 USE grids 00107 USE rotations 00108 USE remap_vars 00109 USE vector 00110 ! 00111 !* ---------------------------- Implicit -------------------------------- 00112 ! 00113 IMPLICIT NONE 00114 ! 00115 !* ---------------------------- Include files --------------------------- 00116 ! 00117 !* ---------------------------- Intent In ------------------------------- 00118 ! 00119 INTEGER(KIND=int_kind), INTENT(in) :: 00120 nlon_src, nlat_src, ! number of longitudes and latitudes on source grid 00121 nlon_dst, nlat_dst, ! number of longitudes and latitudes on target grid 00122 src_size, dst_size, ! number of source/target grid cells 00123 n_srch_bins, ! number of search bins for SCRIP 00124 msk_srcA(src_size), ! grid mask for srcA 00125 msk_srcB(src_size), ! grid mask for srcB 00126 msk_dstA(dst_size), ! grid mask for dstA 00127 msk_dstB(dst_size), ! grid mask for dstB 00128 id_sper, ! number of overlapping points for source grid 00129 id_tper, ! number of overlapping points for target grid 00130 id_scripvoi ! number of neighbour for DISTWGT and GAUSWGT 00131 00132 REAL (KIND=real_kind), INTENT(in) :: 00133 src_arrayI_srcA(src_size), ! field component I on original source grid srcA 00134 src_arrayJ_srcB(src_size), ! field component J on original source grid srcB 00135 lon_srcA(src_size), lat_srcA(src_size), ! lon-/latitudes for source grid srcA 00136 lon_srcB(src_size), lat_srcB(src_size), ! lon-/latitudes for source grid srcB 00137 lon_dstA(dst_size), lat_dstA(dst_size), ! lon-/latitudes for target grid dstA 00138 lon_dstB(dst_size), lat_dstB(dst_size), ! lon-/latitudes for target grid dstB 00139 rl_varmul ! Gaussian variance (for GAUSWGT) 00140 00141 CHARACTER(LEN=8), INTENT(in) :: 00142 grd_name_srcA, grd_name_srcB, ! grid name for the 2 sources grids 00143 grd_name_dstA, grd_name_dstB, ! grid name for the 2 target grids 00144 map_method, ! remapping method 00145 cdgrdtyp, ! source grid type 00146 normalize_opt, ! option for normalization 00147 order, ! order of conservative remapping 00148 rst_type, ! type of scrip search restriction 00149 cd_sper, ! source grid periodicity type 00150 cd_tper ! target grid periodicity type 00151 00152 LOGICAL, INTENT(in) :: 00153 lextrapdone, ! logical, true if EXTRAP done on field 00154 lprojcart ! logical, ture if projection to cartesian should be done 00155 00156 ! 00157 !* ---------------------------- Intent Out ------------------------------- 00158 ! 00159 REAL (KIND=real_kind), INTENT(out):: 00160 dst_arrayI_dstA(dst_size), ! VECTOR_I field on target grid dstA 00161 dst_arrayJ_dstB(dst_size) ! VECTOR_J field on target grid dstB 00162 00163 ! 00164 !* ---------------------------- Local declarations ---------------------- 00165 ! 00166 INTEGER(KIND=int_kind) :: nbr_dst, nbr_comp, ii, jj, icount, kcount, stat, 00167 nc_gridsid, angid, nc_fileid, var_id 00168 INTEGER(KIND=int_kind) :: ilenstr 00169 INTEGER (KIND=int_kind), DIMENSION(dst_size) :: msk_dst 00170 REAL (KIND=real_kind), DIMENSION(src_size) :: src_arrayJ_srcA 00171 REAL (KIND=real_kind), DIMENSION(src_size,3) :: src_array, src_array_buf 00172 REAL (KIND=real_kind), DIMENSION(dst_size) :: lon_dst, lat_dst 00173 REAL (KIND=real_kind), DIMENSION(src_size,2,2) :: mat_loc2sph 00174 REAL (KIND=real_kind), DIMENSION(src_size,3,2) :: mat_sph2car 00175 REAL (KIND=real_kind), DIMENSION(dst_size,3,3) :: mat_car2sph 00176 REAL (KIND=real_kind), DIMENSION(dst_size,2,2) :: mat_sph2loc 00177 REAL (KIND=real_kind), DIMENSION(dst_size,3) :: dst_array, dst_array_buf 00178 00179 CHARACTER(LEN=8) :: grd_name_dst 00180 CHARACTER (char_len) :: cmapping, crmpfile 00181 CHARACTER(LEN=12) :: ang_name 00182 #if defined use_oasis_para || defined use_oasis_cmcc_para 00183 CHARACTER*3 :: cl_indexoa 00184 #endif 00185 00186 LOGICAL :: ll_loc2spher_srcA, ll_spher2loc_dst 00187 00188 ! 00189 !* ---------------------------- Poema verses ---------------------------- 00190 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00191 ! 00192 !*** O. Initialization 00193 ! -------------- 00194 ! -------------- 00195 IF (nlogprt .GE. 2) THEN 00196 WRITE (UNIT = nulou,FMT = *) ' ' 00197 WRITE (UNIT = nulou,FMT = *) & 00198 ' Entering ROUTINE scriprmp_vector - Level 3' 00199 WRITE (UNIT = nulou,FMT = *) & 00200 ' *********************** *******' 00201 WRITE (UNIT = nulou,FMT = *) ' ' 00202 WRITE (UNIT = nulou,FMT = *) 'SCRIP remapping for vector' 00203 WRITE (UNIT = nulou,FMT = *) ' ' 00204 call flush(nulou) 00205 ENDIF 00206 00207 ! 00208 !*** 1. Interpolation towards the same source grid 00209 ! ------------------------------------------ 00210 ! ------------------------------------------ 00211 ! 00212 !* If the source grids of the two vector components are different 00213 !* -> interpolation of field J on grd_srcB towards the same grid as for 00214 !* the field I grd_srcA. 00215 ! 00216 ! 00217 !* After this step the two components are on the same source grid: 00218 !* -> src_arrayI_srcA, src_arrayJ_srcA 00219 ! 00220 00221 IF (grd_name_srcA .NE. grd_name_srcB) THEN 00222 ! Interpolate the field src_arrayJ_srcB from grid srcB towards 00223 ! grid srcA. Resulting field on srcA is src_arrayJ_srcA 00224 CALL scriprmp( & 00225 src_arrayJ_srcA,src_arrayJ_srcB, src_size, src_size, & 00226 msk_srcB, msk_srcA, & 00227 lon_srcB, lat_srcB, nlon_src, nlat_src, & 00228 lon_srcA, lat_srcA, nlon_src, nlat_src, & 00229 map_method, cdgrdtyp, & 00230 id_sper, id_sper, cd_sper, cd_sper, & 00231 grd_name_srcB, grd_name_srcA, & 00232 normalize_opt, order, rst_type, n_srch_bins, & 00233 lextrapdone, rl_varmul, id_scripvoi) 00234 ELSE 00235 ! grid srcA equal grid srcB, no interpolation necessary 00236 src_arrayJ_srcA = src_arrayJ_srcB 00237 ENDIF 00238 00239 ! 00240 !*** 2. Check if the target grids are different 00241 ! If different the interpolation will be done twice 00242 ! ------------------------------------------------- 00243 ! ------------------------------------------------- 00244 00245 IF (grd_name_dstA .EQ. grd_name_dstB) THEN 00246 nbr_dst=1 00247 ELSE 00248 nbr_dst=2 00249 ENDIF 00250 00251 00252 DO ii=1, nbr_dst 00253 00254 IF (ii == 1) THEN 00255 grd_name_dst = grd_name_dstA 00256 lon_dst(:) = lon_dstA(:) 00257 lat_dst(:) = lat_dstA(:) 00258 msk_dst(:) = msk_dstA(:) 00259 ELSE 00260 grd_name_dst = grd_name_dstB 00261 lon_dst(:) = lon_dstB(:) 00262 lat_dst(:) = lat_dstB(:) 00263 msk_dst(:) = msk_dstB(:) 00264 ENDIF 00265 00266 IF (nlogprt .GE. 2) THEN 00267 WRITE (UNIT = nulou,FMT = *) '' 00268 WRITE (UNIT = nulou,FMT = *) & 00269 ' Now making interpolation towards target grid : ', grd_name_dst 00270 CALL FLUSH(nulou) 00271 ENDIF 00272 00273 ! 00274 !*** 3. Calculate and/or read weights, remapping and rotation matrix 00275 ! ------------------------------------------------------------ 00276 ! ------------------------------------------------------------ 00277 ! 00278 !* -- Get the name of the file containing the remapping matrix. 00279 ! 00280 #if defined use_oasis_para || defined use_oasis_cmcc_para 00281 IF (ig_indexoa .le. 9) THEN 00282 WRITE(cl_indexoa,FMT='(I1)') ig_indexoa 00283 ELSE IF (ig_indexoa .le. 99) THEN 00284 WRITE(cl_indexoa,FMT='(I2)') ig_indexoa 00285 ELSE IF (ig_indexoa .le. 999) THEN 00286 WRITE(cl_indexoa,FMT='(I3)') ig_indexoa 00287 ENDIF 00288 #endif 00289 ! 00290 icount = ilenstr(grd_name_srcA,jpeight) 00291 kcount = ilenstr(grd_name_dst, jpeight) 00292 00293 SELECT CASE (map_method) 00294 CASE ('CONSERV') ! conservative remapping 00295 cmapping = & 00296 grd_name_srcA(1:icount)//' to '//grd_name_dst(1:kcount)//' '& 00297 //map_method//' '//normalize_opt(1:4)//' remapping' 00298 #if defined use_oasis_para || defined use_oasis_cmcc_para 00299 IF (ig_indexoa .LE. 9) THEN 00300 crmpfile = & 00301 'rmp_'//grd_name_srcA(1:icount)//'_to_'//grd_name_dst(1:kcount)& 00302 //'_'//map_method(1:7)//'_'//normalize_opt(1:4)// & 00303 '_'//cl_indexoa(1:1)//'.nc' 00304 ELSE IF (ig_indexoa .LE. 99) THEN 00305 crmpfile = & 00306 'rmp_'//grd_name_srcA(1:icount)//'_to_'//grd_name_dst(1:kcount)& 00307 //'_'//map_method(1:7)//'_'//normalize_opt(1:4)// & 00308 '_'//cl_indexoa(1:2)//'.nc' 00309 ELSE IF (ig_indexoa .LE. 999) THEN 00310 crmpfile = & 00311 'rmp_'//grd_name_srcA(1:icount)//'_to_'//grd_name_dst(1:kcount)& 00312 //'_'//map_method(1:7)//'_'//normalize_opt(1:4)// & 00313 '_'//cl_indexoa(1:3)//'.nc' 00314 ENDIF 00315 #else 00316 crmpfile = & 00317 'rmp_'//grd_name_srcA(1:icount)//'_to_'//grd_name_dst(1:kcount)& 00318 //'_'//map_method(1:7)//'_'//normalize_opt(1:4)//'.nc' 00319 #endif 00320 ! 00321 CASE DEFAULT 00322 cmapping = & 00323 grd_name_srcA(1:icount)//' to '//grd_name_dst(1:kcount)//' '& 00324 //map_method//' remapping' 00325 #if defined use_oasis_para || defined use_oasis_cmcc_para 00326 IF (ig_indexoa .LE. 9) THEN 00327 crmpfile = & 00328 'rmp_'//grd_name_srcA(1:icount)//'_to_'//grd_name_dst(1:kcount)& 00329 //'_'//map_method(1:7)//'_'//cl_indexoa(1:1)//'.nc' 00330 ELSE IF (ig_indexoa .le. 99) THEN 00331 crmpfile = & 00332 'rmp_'//grd_name_srcA(1:icount)//'_to_'//grd_name_dst(1:kcount)& 00333 //'_'//map_method(1:7)//'_'//cl_indexoa(1:2)//'.nc' 00334 ELSE IF (ig_indexoa .le. 999) THEN 00335 crmpfile = & 00336 'rmp_'//grd_name_srcA(1:icount)//'_to_'//grd_name_dst(1:kcount)& 00337 //'_'//map_method(1:7)//'_'//cl_indexoa(1:3)//'.nc' 00338 ENDIF 00339 #else 00340 crmpfile = & 00341 'rmp_'//grd_name_srcA(1:icount)//'_to_'//grd_name_dst(1:kcount)& 00342 //'_'//map_method(1:7)//'.nc' 00343 #endif 00344 END SELECT 00345 00346 IF (nlogprt .GE. 2) THEN 00347 WRITE (UNIT = nulou,FMT = *) & 00348 ' SCRIP filename : ', crmpfile 00349 CALL FLUSH(nulou) 00350 ENDIF 00351 00352 ! 00353 !*** 3.1 Create and/or open the file with remapping matrix 00354 ! ------------------------------------------------- 00355 ! 00356 stat = NF_OPEN(crmpfile, NF_WRITE, nc_fileid) 00357 IF (stat .NE. NF_NOERR) THEN 00358 ! 00359 !* -- The file does not exist. The remapping matrix has to be created 00360 ! for remapping from source grd_name_srcA to target grd_name_dst 00361 ! 00362 CALL calc_remap_matrix( & 00363 msk_srcA(:), msk_dst(:), & 00364 src_size, dst_size, & 00365 lon_srcA(:), lat_srcA(:), lon_dst(:), lat_dst(:), & 00366 nlon_src, nlat_src, nlon_dst, nlat_dst, & 00367 grd_name_srcA, grd_name_dst, & 00368 map_method, normalize_opt, cdgrdtyp, & 00369 id_sper, cd_sper, & 00370 id_tper, cd_tper, & 00371 rst_type, & 00372 n_srch_bins, crmpfile, & 00373 cmapping, lextrapdone, & 00374 rl_varmul, id_scripvoi ) 00375 00376 stat = NF_OPEN(crmpfile, NF_WRITE, nc_fileid) 00377 00378 ELSE 00379 00380 !* -- The file with remapping matrix existe already 00381 ! 00382 IF (nlogprt .GE. 2) THEN 00383 WRITE (UNIT = nulou,FMT = *) & 00384 'SCRIP file exists with remapping matrix src to dst grid' 00385 WRITE (UNIT = nulou,FMT = *) ' ' 00386 CALL FLUSH(nulou) 00387 ENDIF 00388 00389 END IF !IF (stat .NE. NF_NOERR) 00390 00391 00392 !*** 3.2. If necessary calculate matrix for rotation on srcA grid local to sperical 00393 ! -------------------------------------------------------------------------- 00394 ! 00395 !* Get angles from file grids.nc. 00396 ! If they exist a rotation from local to spherical referntial will be done, 00397 ! if no angles the grid is considered as a spheric referential 00398 ! 00399 CALL hdlerr(NF_OPEN & 00400 ('grids.nc',NF_NOWRITE, nc_gridsid), 'scriprmp_vector') 00401 ang_name = grd_name_srcA(1:icount)//'.ang' 00402 stat = NF_INQ_VARID(nc_gridsid, ang_name, var_id) 00403 00404 CALL hdlerr(NF_CLOSE(nc_gridsid), 'scriprmp_vector') 00405 00406 IF(stat == NF_NOERR) THEN 00407 ll_loc2spher_srcA = .true. 00408 ! 00409 !* Rotaion from local to spherical referential on grd_srcA is needed 00410 !* Check if the rotation matrix Mat_src_loc2spher exists in the remapping file 00411 ! 00412 stat = NF_INQ_VARID(nc_fileid, 'Mat_src_loc2sph', var_id) 00413 00414 IF (stat == NF_NOERR) THEN 00415 !* Read the matrix already existing 00416 IF (ll_single) THEN 00417 CALL hdlerr(NF_GET_VAR_REAL & 00418 (nc_fileid, var_id, mat_loc2sph(:,:,:)), 'scriprmp_vector') 00419 ELSE 00420 CALL hdlerr(NF_GET_VAR_DOUBLE & 00421 (nc_fileid, var_id, mat_loc2sph(:,:,:)), 'scriprmp_vector') 00422 ENDIF 00423 ELSE 00424 !* The matrix has to be calculated and written to the file 00425 !* CALL loc2spher() writes the matrix to netcdf file nc_fileid 00426 ! and returns the matrix mat_loc2sph 00427 CALL loc2spher(grd_name_srcA, src_size, mat_loc2sph(:,:,:), nc_fileid) 00428 ENDIF 00429 00430 ELSE ! The angles does not exist in grids.nc 00431 ll_loc2spher_srcA = .false. 00432 WRITE (UNIT = nulou,FMT = *) & 00433 ' WARNING' 00434 WRITE (UNIT = nulou,FMT = *) & 00435 '***************' 00436 WRITE (UNIT = nulou,FMT = *) & 00437 'file grids.nc contains no angle information' 00438 WRITE (UNIT = nulou,FMT = *) & 00439 'source grid '// grd_name_srcA// & 00440 ' will be considered as a spheric referential' 00441 CALL FLUSH(nulou) 00442 END IF 00443 ! 00444 !*** 3.3 Calculate matrix for projection to cartesian if demanded by user option PROJCART 00445 ! 00446 ! -> calculate matrix from spherical to cartesian referntials on source grid srcA 00447 ! -> calculate matrix from cartesian to spherical referntials on target grid dst 00448 ! 00449 IF ( lprojcart ) THEN 00450 00451 !* On source grid 00452 00453 stat = NF_INQ_VARID(nc_fileid, 'Mat_src_sph2car', var_id) 00454 IF (stat == NF_NOERR) THEN 00455 !* Read the matrix already existing 00456 IF (ll_single) THEN 00457 CALL hdlerr(NF_GET_VAR_REAL & 00458 (nc_fileid, var_id, mat_sph2car(:,:,:)), 'scriprmp_vector') 00459 ELSE 00460 CALL hdlerr(NF_GET_VAR_DOUBLE & 00461 (nc_fileid, var_id, mat_sph2car(:,:,:)), 'scriprmp_vector') 00462 ENDIF 00463 ELSE 00464 !* Calculate the matrix mat_sph2car and write to file nc_fileid 00465 ! 00466 CALL spher2car(lon_srcA(:), lat_srcA(:), src_size, mat_sph2car(:,:,:), nc_fileid) 00467 00468 ENDIF 00469 00470 !* On target grid 00471 00472 stat = NF_INQ_VARID(nc_fileid, 'Mat_dst_car2sph', var_id) 00473 IF (stat == NF_NOERR) THEN 00474 !* Read the matrix already existing 00475 IF (ll_single) THEN 00476 CALL hdlerr(NF_GET_VAR_REAL & 00477 (nc_fileid, var_id, mat_car2sph(:,:,:)), 'scriprmp_vector') 00478 ELSE 00479 CALL hdlerr(NF_GET_VAR_DOUBLE & 00480 (nc_fileid, var_id, mat_car2sph(:,:,:)), 'scriprmp_vector') 00481 ENDIF 00482 ELSE 00483 !* Calculate the matrix mat_car2sph and write to file nc_fileid 00484 ! 00485 CALL car2spher(lon_dst(:), lat_dst(:), dst_size, mat_car2sph(:,:,:), nc_fileid) 00486 00487 ENDIF 00488 00489 ENDIF !IF ( lprojcart ) 00490 00491 ! 00492 !*** 3.4 If necessary calculate matrix for rotation on target grid sperical to local 00493 ! --------------------------------------------------------------------------- 00494 ! 00495 !* Get angles from file grids.nc. 00496 ! If they exist a rotation from spherical to local referntial will be done, 00497 ! if no angles the grid is considered as a spheric referential 00498 ! 00499 CALL hdlerr(NF_OPEN & 00500 ('grids.nc',NF_NOWRITE, nc_gridsid), 'scriprmp_vector') 00501 ang_name = grd_name_dst(1:kcount)//'.ang' 00502 stat = NF_INQ_VARID(nc_gridsid, ang_name, angid) 00503 00504 CALL hdlerr(NF_CLOSE(nc_gridsid), 'scriprmp_vector') 00505 00506 IF (stat == NF_NOERR) THEN 00507 ll_spher2loc_dst = .true. 00508 ! 00509 !* Rotaion from spherical to local referential on grd_dst is needed 00510 ! Check if the rotation matrix Mat_dst_sph2loc exists in the remapping file 00511 ! 00512 stat = NF_INQ_VARID(nc_fileid, 'Mat_dst_sph2loc', var_id) 00513 00514 IF (stat == NF_NOERR) THEN 00515 !* Read the matrix already existing 00516 IF (ll_single) THEN 00517 CALL hdlerr(NF_GET_VAR_REAL & 00518 (nc_fileid, var_id, mat_sph2loc(:,:,:)), 'scriprmp_vector') 00519 ELSE 00520 CALL hdlerr(NF_GET_VAR_DOUBLE & 00521 (nc_fileid, var_id, mat_sph2loc(:,:,:)), 'scriprmp_vector') 00522 ENDIF 00523 ELSE 00524 !* The matrix has to be calculated and written to the file 00525 !* CALL spher2loc() writes the matrix to netcdf file nc_fileid 00526 ! and returns the matrix mat_sph2loc 00527 CALL spher2loc(grd_name_dst, dst_size, mat_sph2loc(:,:,:), nc_fileid) 00528 00529 ENDIF 00530 00531 ELSE ! The angles does not exist in grids.nc 00532 ll_spher2loc_dst = .false. 00533 WRITE (UNIT = nulou,FMT = *) & 00534 ' WARNING' 00535 WRITE (UNIT = nulou,FMT = *) & 00536 '***************' 00537 WRITE (UNIT = nulou,FMT = *) & 00538 'file grids.nc contains no angle information' 00539 WRITE (UNIT = nulou,FMT = *) & 00540 'target grid '// grd_name_dst// & 00541 ' will be considered as a spheric referential' 00542 CALL FLUSH(nulou) 00543 00544 END IF ! IF (stat == NF_NOERR) 00545 00546 ! 00547 !* End of calculation of weight and transformation matrix 00548 !* All matrix do now exist 00549 ! 00550 00551 ! 00552 !*** 4. Rotation and interpolation of the two components on grd_srcA towards 00553 ! the target grid grd_dst 00554 ! ------------------------------------------------------------------- 00555 ! ------------------------------------------------------------------- 00556 ! 00557 !* save the both source components in src_array(:,:) 00558 src_array(:,1) = src_arrayI_srcA(:) 00559 src_array(:,2) = src_arrayJ_srcA(:) 00560 src_array(:,3) = 0 00561 nbr_comp=2 00562 ! 00563 !*** 4.1 Rotate the the vector field components on srcA to spherical referential 00564 ! 00565 IF ( ll_loc2spher_srcA ) THEN 00566 00567 src_array(:,1) = mat_loc2sph(:,1,1) * src_arrayI_srcA(:) + & 00568 mat_loc2sph(:,1,2) * src_arrayJ_srcA(:) 00569 00570 src_array(:,2) = mat_loc2sph(:,2,1) * src_arrayI_srcA(:) + & 00571 mat_loc2sph(:,2,2) * src_arrayJ_srcA(:) 00572 00573 ENDIF 00574 #ifdef __DEBUG 00575 ! The components in spherical referential are saved to vector_debugXXXXXX.nc file. 00576 CALL write_src_array_spheric(src_array(:,1:2), src_size, nlon_src, nlat_src, & 00577 grd_name_srcA, grd_name_dst) 00578 #endif 00579 00580 ! 00581 !*** 4.2 If PROJCART, rotate the the vector field components to cartesian referential 00582 ! 00583 IF (lprojcart) THEN 00584 00585 nbr_comp=3 00586 00587 src_array_buf(:,1) = & 00588 mat_sph2car(:,1,1) * src_array(:,1) + & 00589 mat_sph2car(:,1,2) * src_array(:,2) 00590 00591 src_array_buf(:,2) = & 00592 mat_sph2car(:,2,1) * src_array(:,1) + & 00593 mat_sph2car(:,2,2) * src_array(:,2) 00594 00595 src_array_buf(:,3) = & 00596 mat_sph2car(:,3,1) * src_array(:,1) + & 00597 mat_sph2car(:,3,2) * src_array(:,2) 00598 00599 src_array(:,1)= src_array_buf(:,1) 00600 src_array(:,2)= src_array_buf(:,2) 00601 src_array(:,3)= src_array_buf(:,3) 00602 00603 ENDIF 00604 ! 00605 !*** 4.3 Interpolate scalar the 2 or 3 compontents src_array(:,:) 00606 ! from the source grd_srcA to target grd_dst, result in dst_array(:,:) 00607 ! -------------------------------------------------------------------- 00608 00609 CALL remap_vector_comps (dst_array(:,1:nbr_comp), src_array(:,1:nbr_comp), & 00610 nbr_comp, nc_fileid, & 00611 map_method, cdgrdtyp, order, & 00612 msk_srcA(:), & 00613 id_sper, cd_sper, & 00614 nlon_src, nlat_src, lon_srcA(:), lat_srcA(:), src_size, dst_size) 00615 00616 ! 00617 !** Close remapping file 00618 ! 00619 CALL hdlerr(NF_CLOSE(nc_fileid), 'scriprmp') 00620 00621 ! 00622 !*** 4.4 If PROJCART, rotate dst_array from cartesian to spherical referential 00623 ! --------------------------------------------------------------------- 00624 00625 IF (lprojcart) THEN 00626 00627 #ifdef __DEBUG 00628 ! The components in cartesian referential are saved to vector_debugXXXXXX.nc file 00629 CALL write_cartesian_components(src_array(:,1:3), dst_array(:,1:3), src_size, & 00630 dst_size, nlon_src, nlat_src, nlon_dst, nlat_dst, grd_name_srcA, & 00631 grd_name_dst) 00632 #endif 00633 00634 dst_array_buf(:,1) = & 00635 mat_car2sph(:,1,1) * dst_array(:,1) + & 00636 mat_car2sph(:,1,2) * dst_array(:,2) + & 00637 mat_car2sph(:,1,3) * dst_array(:,3) 00638 00639 dst_array_buf(:,2) = & 00640 mat_car2sph(:,2,1) * dst_array(:,1) + & 00641 mat_car2sph(:,2,2) * dst_array(:,2) + & 00642 mat_car2sph(:,2,3) * dst_array(:,3) 00643 00644 dst_array_buf(:,3) = & 00645 mat_car2sph(:,3,1) * dst_array(:,1) + & 00646 mat_car2sph(:,3,2) * dst_array(:,2) + & 00647 mat_car2sph(:,3,3) * dst_array(:,3) 00648 00649 dst_array(:,1) = dst_array_buf(:,1) 00650 dst_array(:,2) = dst_array_buf(:,2) 00651 dst_array(:,3) = dst_array_buf(:,3) 00652 00653 ENDIF 00654 00655 ! 00656 !*** 4.5 If the target grid is not in a spherical referential 00657 ! then rotate the spherical components towords locals 00658 ! ---------------------------------------------------- 00659 00660 IF (ll_spher2loc_dst) THEN 00661 00662 #ifdef __DEBUG 00663 ! The components in spherical referential are saved to vector_debugXXXXXX.nc file 00664 ! The field K dst_array(:,3) should be close to zero 00665 CALL write_dst_array_spheric(dst_array(:,1:2), dst_size, nlon_dst, nlat_dst, & 00666 grd_name_srcA, grd_name_dst) 00667 #endif 00668 00669 dst_array_buf(:,1) = & 00670 mat_sph2loc(:,1,1) * dst_array(:,1) + & 00671 mat_sph2loc(:,1,2) * dst_array(:,2) 00672 00673 dst_array_buf(:,2) = & 00674 mat_sph2loc(:,2,1) * dst_array(:,1) + & 00675 mat_sph2loc(:,2,2) * dst_array(:,2) 00676 00677 dst_array(:,1) = dst_array_buf(:,1) 00678 dst_array(:,2) = dst_array_buf(:,2) 00679 00680 ENDIF 00681 00682 !* The dst_array is now the I and J components in the local referential 00683 !* on target grid grd_dst 00684 ! 00685 00686 00687 ! 00688 !*** 5. Save the the components that correspond to VECTOR_I and VECTOR_J grid 00689 ! --------------------------------------------------------------------- 00690 ! --------------------------------------------------------------------- 00691 00692 IF (ii == 1 .AND. nbr_dst == 1) THEN 00693 ! There is only one target grid. Save both components calculated above. 00694 dst_arrayI_dstA(:) = dst_array(:,1) 00695 dst_arrayJ_dstB(:) = dst_array(:,2) 00696 ELSEIF (ii == 1) THEN 00697 ! There are 2 target grid and the one just calculated is the grid for VECTOR_I 00698 dst_arrayI_dstA(:) = dst_array(:,1) 00699 ELSE 00700 ! There are 2 target grids and the one just calculated is the grid for VECTOR_J 00701 dst_arrayJ_dstB(:) = dst_array(:,2) 00702 ENDIF 00703 00704 IF (nlogprt .GE. 2) THEN 00705 write(nulou,*)' src_size = ',src_size,'dst_size = ',dst_size 00706 write(nulou,*)' Now finish interpolation towards target grid : ', grd_name_dst 00707 call flush(nulou) 00708 ENDIF 00709 00710 END DO 00711 00712 ! 00713 !*** 6. Test if there are points on the latitudes 90deg north and 90deg south. 00714 ! If there are points on these latitudes, the averege for these points 00715 ! are calculated and distributed. 00716 ! ----------------------------------------------------------------------- 00717 ! ----------------------------------------------------------------------- 00718 00719 CALL check_points_at_poles (dst_arrayI_dstA(:), dst_arrayJ_dstB(:), lat_dstA(:), & 00720 lat_dstB(:), grd_name_dstA, grd_name_dstB, dst_size) 00721 00722 ! 00723 !* End of routine 00724 ! 00725 IF (nlogprt .GE. 2) THEN 00726 WRITE (UNIT = nulou,FMT = *) ' ' 00727 WRITE (UNIT = nulou,FMT = *) & 00728 ' --------- End of routine scriprmp_vector ---------' 00729 CALL FLUSH (nulou) 00730 ENDIF 00731 00732 END SUBROUTINE scriprmp_vector
1.7.4 
