|
Oasis3 4.0.2
|
|
Oasis3 4.0.2
|
00001 MODULE vector 00002 !----------------------------------------------------------------------- 00003 ! 00004 !** Description: 00005 ! ------------ 00006 ! This module contains subroutines needed for scriprmp_vector 00007 ! 00008 !** History: 00009 ! -------- 00010 ! Created 02/2006 J. Ghattas 00011 ! 00012 ! 00013 USE kinds_mod 00014 USE constants 00015 USE mod_parameter 00016 USE grids 00017 USE mod_unitncdf 00018 00019 IMPLICIT NONE 00020 00021 #include <netcdf.inc> 00022 00023 00024 CONTAINS 00025 00026 ! 00027 !** ----------------------------------------------------------------------- 00028 !** ----------------------------------------------------------------------- 00029 ! 00030 00031 SUBROUTINE calc_remap_matrix( & 00032 src_mask, dst_mask, & 00033 src_size, dst_size, & 00034 lon_srcA, lat_srcA, lon_dst, lat_dst, & 00035 nlon_src, nlat_src, nlon_dst, nlat_dst, & 00036 grd_name_srcA, grd_name_dst, & 00037 map_method, n_opt, cdgrdtyp, & 00038 id_sper, cd_sper, & 00039 id_tper, cd_tper, & 00040 rst_type, & 00041 n_srch_bins, crmpfile, & 00042 cmapping, lextrapdone, & 00043 rl_varmul, id_scripvoi ) 00044 00045 ! 00046 !** Description : 00047 ! ------------- 00048 ! Calcutation of remapping matrix for interpolation between qource grid grd_name_src 00049 ! and target grid grd_name_dst. Remapping matrix written to file crmpfile which must 00050 ! not exist when the this subroutine is called. 00051 ! 00052 IMPLICIT NONE 00053 ! 00054 !** Input variables 00055 ! --------------- 00056 INTEGER(KIND=int_kind), INTENT(in) :: src_size, dst_size, 00057 nlon_src, nlat_src, nlon_dst, nlat_dst, 00058 id_sper, id_tper, n_srch_bins, id_scripvoi 00059 INTEGER(KIND=int_kind), DIMENSION(src_size), INTENT(in) :: src_mask 00060 INTEGER(KIND=int_kind), DIMENSION(dst_size), INTENT(in) :: dst_mask 00061 REAL(KIND=real_kind), DIMENSION(src_size), INTENT(in) :: lon_srcA, lat_srcA 00062 REAL(KIND=real_kind), DIMENSION(dst_size), INTENT(in) :: lon_dst, lat_dst 00063 REAL(KIND=real_kind), INTENT(in) :: rl_varmul 00064 CHARACTER(LEN=8), INTENT(in) :: grd_name_srcA, grd_name_dst, n_opt, rst_type 00065 CHARACTER(LEN=8), INTENT(in) :: map_method, cdgrdtyp, cd_sper, cd_tper 00066 CHARACTER(char_len), INTENT(in) :: crmpfile, cmapping 00067 LOGICAL, INTENT(in) :: lextrapdone 00068 00069 ! 00070 !** Local variables 00071 ! --------------- 00072 INTEGER(KIND=int_kind), DIMENSION(src_size) :: sou_mask 00073 INTEGER(KIND=int_kind), DIMENSION(dst_size) :: tgt_mask 00074 CHARACTER(LEN=8) :: normalize_opt 00075 LOGICAL :: lfracnnei 00076 INTEGER(KIND=int_kind) :: ncrn, src_rank, dst_rank, src_dims(2), dst_dims(2) 00077 REAL (KIND=real_kind), DIMENSION(:,:), ALLOCATABLE :: src_corner_lon, src_corner_lat, 00078 dst_corner_lon, dst_corner_lat 00079 ! 00080 CHARACTER(LEN=8) :: cl_tgt ! indicates target grid 00081 ! 00082 !** ---------------------------------------------------------------------------- 00083 ! 00084 ! 00085 IF (nlogprt .GE. 2) THEN 00086 WRITE (UNIT = nulou,FMT = *)' ' 00087 WRITE (UNIT = nulou,FMT = *)'Entering routine calc_remap_matrix' 00088 CALL FLUSH(nulou) 00089 ENDIF 00090 ! 00091 ! 00092 !** -- Setting of the mask for the source and the target grid 00093 ! 00094 00095 sou_mask(:) = 0 00096 tgt_mask(:) = 0 00097 WHERE (src_mask .eq. 1) 00098 sou_mask = 0 00099 END WHERE 00100 WHERE (src_mask .eq. 0) 00101 sou_mask = 1 00102 END WHERE 00103 00104 WHERE (dst_mask .eq. 1) 00105 tgt_mask = 0 00106 END WHERE 00107 WHERE (dst_mask .eq. 0) 00108 tgt_mask = 1 00109 END WHERE 00110 00111 ! 00112 !** -- Calculate SCRIP remapping matrix : source to destination grid 00113 ! 00114 normalize_opt = n_opt 00115 SELECT CASE (normalize_opt) 00116 CASE ('FRACNNEI') 00117 normalize_opt = 'FRACAREA' 00118 lfracnnei = .true. 00119 END SELECT 00120 00121 IF (nlogprt .GE. 2) THEN 00122 WRITE (UNIT = nulou,FMT = *) & 00123 ' Calculation of SCRIP remapping matrix: method = ', & 00124 map_method 00125 WRITE (UNIT = nulou,FMT = *) ' ' 00126 call flush(nulou) 00127 ENDIF 00128 ! 00129 !** -- Get grid cell corners for conservative remapping 00130 ! 00131 ncrn = 4. 00132 ALLOCATE(src_corner_lon(ncrn,src_size), & 00133 src_corner_lat(ncrn,src_size),& 00134 dst_corner_lon(ncrn,dst_size),& 00135 dst_corner_lat(ncrn,dst_size)) 00136 src_corner_lon(:,:)= 0.0 00137 src_corner_lat(:,:)= 0.0 00138 dst_corner_lon(:,:)= 0.0 00139 dst_corner_lat(:,:)= 0.0 00140 00141 IF (map_method == 'CONSERV') THEN 00142 CALL corners(nlon_src, nlat_src, ncrn, & 00143 lon_srcA, lat_srcA, & 00144 grd_name_srcA, cdgrdtyp, id_sper, cd_sper, & 00145 src_corner_lon, src_corner_lat) 00146 cl_tgt='TARGETGR' 00147 CALL corners(nlon_dst, nlat_dst, ncrn, & 00148 lon_dst, lat_dst, & 00149 grd_name_dst, cl_tgt, id_tper, cd_tper, & 00150 dst_corner_lon, dst_corner_lat) 00151 ENDIF 00152 ! 00153 !** -- Initialization of grid arrays for SCRIP 00154 ! 00155 src_dims(1) = nlon_src 00156 src_dims(2) = nlat_src 00157 dst_dims(1) = nlon_dst 00158 dst_dims(2) = nlat_dst 00159 src_rank = 2 00160 dst_rank = 2 00161 00162 CALL grid_init(map_method, rst_type, n_srch_bins,& 00163 src_size, dst_size, src_dims(:), dst_dims(:),& 00164 src_rank, dst_rank, ncrn, ncrn,& 00165 sou_mask(:), tgt_mask(:), grd_name_srcA, grd_name_dst,& 00166 lat_srcA(:), lon_srcA(:), lat_dst(:), lon_dst(:),& 00167 src_corner_lat(:,:), src_corner_lon(:,:),& 00168 dst_corner_lat(:,:), dst_corner_lon(:,:)) 00169 ! 00170 !** -- Calculation of weights and addresses using SCRIP-library 00171 ! 00172 CALL scrip(crmpfile, cmapping, map_method, normalize_opt, & 00173 lextrapdone, rl_varmul, id_scripvoi) 00174 00175 IF (map_method == 'CONSERV') THEN 00176 DEALLOCATE(src_corner_lon, src_corner_lat,& 00177 dst_corner_lon, dst_corner_lat) 00178 ENDIF 00179 00180 ! 00181 IF (nlogprt .GE. 2) THEN 00182 WRITE (UNIT = nulou,FMT = *)' ' 00183 WRITE (UNIT = nulou,FMT = *)'Leaving routine calc_remap_matrix' 00184 CALL FLUSH(nulou) 00185 ENDIF 00186 ! 00187 END SUBROUTINE calc_remap_matrix 00188 00189 ! 00190 !** ----------------------------------------------------------------------- 00191 !** ----------------------------------------------------------------------- 00192 ! 00193 00194 SUBROUTINE remap_vector_comps (dst_array, src_array, & 00195 nbr_comp, nc_fileid, & 00196 map_method, cdgrdtyp, order, & 00197 src_mask, & 00198 id_sper, cd_sper, & 00199 nlon_src, nlat_src, src_lon, src_lat, src_size, dst_size) 00200 00201 ! 00202 !** Description : 00203 ! ------------- 00204 ! Interpolation of 2 or 3 components in src_array(:,:) on source grid using 00205 ! remapping matrix allready existing in remappingile. Resulting fields in 00206 ! dst_array on target grid grd_name_dst. Remapping matrix is read from file 00207 ! already opend with id nc_fileid. 00208 ! 00209 IMPLICIT NONE 00210 ! 00211 !** Input variables 00212 ! --------------- 00213 INTEGER (KIND=int_kind), INTENT(in) :: nbr_comp, nc_fileid, id_sper, nlon_src, 00214 nlat_src, src_size, dst_size 00215 REAL (KIND=real_kind), DIMENSION (src_size, nbr_comp), INTENT(in) :: src_array 00216 REAL (KIND=real_kind), DIMENSION (src_size), INTENT(in) :: src_lon, src_lat 00217 CHARACTER (LEN=8), INTENT(in) :: map_method, cdgrdtyp, order, cd_sper 00218 INTEGER (KIND=int_kind), DIMENSION (src_size), INTENT(in) :: src_mask 00219 ! 00220 !** Output variables 00221 ! ---------------- 00222 REAL (KIND=real_kind), DIMENSION (dst_size, nbr_comp), INTENT(out) :: dst_array 00223 ! 00224 !** Local variables 00225 ! --------------- 00226 INTEGER (KIND=int_kind) :: num_wgts, dimid1, num_links, ib, n, varid 00227 00228 INTEGER (KIND=int_kind), DIMENSION (:), ALLOCATABLE :: src_addr, dst_addr 00229 INTEGER (KIND=int_kind), DIMENSION (src_size) :: sou_mask 00230 00231 REAL (KIND=real_kind), DIMENSION (:,:), ALLOCATABLE :: weights 00232 00233 REAL (KIND=real_kind), DIMENSION (:), ALLOCATABLE :: gradient_lat, gradient_lon, 00234 gradient_i, gradient_j, gradient_ij 00235 00236 REAL (KIND=real_kind), DIMENSION (dst_size) :: dst_area, dst_frac 00237 00238 CHARACTER(LEN=11) :: 00239 csrcadd, ! string for source grid addresses 00240 cdstadd ! string for destination grid addresses 00241 CHARACTER(LEN=13) :: 00242 cdstare, ! string for destination grid area 00243 cdstfra ! string for destination grid frac 00244 CHARACTER(LEN=12) :: 00245 cweight ! string for weights 00246 REAL (KIND=real_kind), DIMENSION (dst_size, nbr_comp) :: weightot 00247 LOGICAL :: ll_weightot 00248 00249 ! 00250 !** ---------------------------------------------------------------------------- 00251 ! 00252 IF (nlogprt .GE. 2) THEN 00253 WRITE (UNIT = nulou,FMT = *)' ' 00254 WRITE (UNIT = nulou,FMT = *)'Entering routine remap_vector_comps' 00255 CALL FLUSH(nulou) 00256 ENDIF 00257 ! 00258 !* Read weights and addresses 00259 ! 00260 !** Get number of weights 00261 ! 00262 SELECT CASE (map_method) 00263 CASE ('CONSERV') ! conservative remapping 00264 num_wgts = 3. 00265 CASE ('BILINEAR') ! bilinear remapping 00266 num_wgts = 1. 00267 CASE ('BICUBIC') ! bicubic remapping 00268 IF (cdgrdtyp .eq. 'LR') THEN ! logically rectangular 00269 num_wgts = 4. 00270 ELSE ! reduced 00271 num_wgts=1. 00272 ENDIF 00273 CASE ('DISTWGT') ! distance weighted averaging 00274 num_wgts = 1. 00275 CASE ('GAUSWGT') ! distance gaussian weighted averaging 00276 num_wgts = 1. 00277 END SELECT 00278 00279 ! 00280 !** Setting of the mask for the source grid 00281 ! 00282 sou_mask(:) = 0 00283 00284 WHERE (src_mask .eq. 1) 00285 sou_mask = 0 00286 END WHERE 00287 WHERE (src_mask .eq. 0) 00288 sou_mask = 1 00289 END WHERE 00290 00291 ! 00292 !** Character strings of weights and addresses 00293 ! 00294 csrcadd = 'src_address' 00295 cdstadd = 'dst_address' 00296 cweight = 'remap_matrix' 00297 cdstare = 'dst_grid_area' 00298 cdstfra = 'dst_grid_frac' 00299 ! 00300 !** Get matrix size 00301 ! 00302 CALL hdlerr(NF_INQ_DIMID & 00303 (nc_fileid, 'num_links', dimid1), 'scriprmp_vector') 00304 CALL hdlerr(NF_INQ_DIMLEN & 00305 (nc_fileid, dimid1, num_links), 'scriprmp_vector') 00306 ! 00307 !** Array allocation 00308 ! 00309 ALLOCATE (src_addr(num_links), dst_addr(num_links), & 00310 weights(num_wgts,num_links)) 00311 ! 00312 !** Read source grid addresses and weights 00313 ! 00314 CALL hdlerr(NF_INQ_VARID & 00315 (nc_fileid, csrcadd, varid), 'scriprmp_vector') 00316 CALL hdlerr(NF_GET_VAR_INT & 00317 (nc_fileid, varid, src_addr), 'scriprmp_vector') 00318 CALL hdlerr(NF_INQ_VARID & 00319 (nc_fileid, cdstadd, varid), 'scriprmp_vector') 00320 CALL hdlerr(NF_GET_VAR_INT & 00321 (nc_fileid, varid, dst_addr), 'scriprmp_vector') 00322 CALL hdlerr(NF_INQ_VARID & 00323 (nc_fileid, cweight, varid), 'scriprmp_vector') 00324 IF (ll_single) THEN 00325 CALL hdlerr(NF_GET_VAR_REAL & 00326 (nc_fileid, varid, weights), 'scriprmp_vector') 00327 ELSE 00328 CALL hdlerr(NF_GET_VAR_DOUBLE & 00329 (nc_fileid, varid, weights), 'scriprmp_vector') 00330 ENDIF 00331 CALL hdlerr(NF_INQ_VARID & 00332 (nc_fileid, cdstare, varid), 'scriprmp_vector') 00333 IF (ll_single) THEN 00334 CALL hdlerr(NF_GET_VAR_REAL & 00335 (nc_fileid, varid, dst_area), 'scriprmp_vector') 00336 ELSE 00337 CALL hdlerr(NF_GET_VAR_DOUBLE & 00338 (nc_fileid, varid, dst_area), 'scriprmp_vector') 00339 ENDIF 00340 CALL hdlerr(NF_INQ_VARID & 00341 (nc_fileid, cdstfra, varid), 'scriprmp_vector') 00342 IF (ll_single) THEN 00343 CALL hdlerr(NF_GET_VAR_REAL & 00344 (nc_fileid, varid, dst_frac), 'scriprmp_vector') 00345 ELSE 00346 CALL hdlerr(NF_GET_VAR_DOUBLE & 00347 (nc_fileid, varid, dst_frac), 'scriprmp_vector') 00348 ENDIF 00349 00350 ! 00351 !** Do the matrix multiplication for the 2 or 3 components 00352 ! 00353 ll_weightot = .false. 00354 weightot(:,:) = 0.0 00355 dst_array(:,:) = 0. 00356 00357 SELECT CASE (map_method) 00358 00359 CASE ('CONSERV') ! conservative remapping 00360 SELECT CASE (order) 00361 CASE ('FIRST') ! first order remapping 00362 DO ib=1,nbr_comp 00363 DO n = 1, num_links 00364 IF (src_addr(n) .NE. 0) THEN 00365 dst_array(dst_addr(n),ib) = dst_array(dst_addr(n),ib)& 00366 + weights(1,n) * src_array(src_addr(n),ib) 00367 weightot(dst_addr(n),ib) = weightot(dst_addr(n),ib) & 00368 + weights(1,n) 00369 ENDIF 00370 END DO 00371 END DO 00372 CASE ('SECOND') ! second order remapping (including gradients) 00373 IF (cdgrdtyp .ne. 'LR') THEN 00374 WRITE (UNIT = nulou,FMT = *) & 00375 'Field gradient cannot be calculated' 00376 WRITE (UNIT = nulou,FMT = *) & 00377 'by Oasis as grid is not logically rectangular' 00378 CALL HALTE('STOP in scriprmp (CONSERV)') 00379 ENDIF 00380 ALLOCATE(gradient_lat(src_size), gradient_lon(src_size)) 00381 00382 DO ib=1, nbr_comp 00383 CALL gradient(nlon_src, nlat_src, src_array(:,ib), & 00384 sou_mask, src_lat, src_lon, id_sper, cd_sper, & 00385 gradient_lat, gradient_lon) 00386 00387 DO n = 1, num_links 00388 IF (src_addr(n) .NE. 0) THEN 00389 dst_array(dst_addr(n),ib) = dst_array(dst_addr(n),ib)& 00390 + weights(1,n) * src_array(src_addr(n),ib) & 00391 + weights(2,n) * gradient_lat(src_addr(n)) & 00392 + weights(3,n) * gradient_lon(src_addr(n)) 00393 weightot(dst_addr(n),ib) = weightot(dst_addr(n),ib) & 00394 + weights(1,n) + weights(2,n) + weights(3,n) 00395 ENDIF 00396 ENDDO 00397 END DO 00398 DEALLOCATE(gradient_lat, gradient_lon) 00399 END SELECT! order 00400 CASE ('BILINEAR') ! bilinear remapping 00401 DO ib=1, nbr_comp 00402 DO n = 1, num_links 00403 IF (src_addr(n) .NE. 0) THEN 00404 dst_array(dst_addr(n),ib) = dst_array(dst_addr(n),ib) & 00405 + weights(1,n) * src_array(src_addr(n),ib) 00406 weightot(dst_addr(n),ib) = weightot(dst_addr(n),ib) & 00407 + weights(1,n) 00408 ENDIF 00409 ENDDO 00410 END DO 00411 CASE ('BICUBIC') ! bicubic remapping 00412 SELECT CASE (cdgrdtyp) ! 00413 CASE ('LR') ! logically rectangular 00414 ALLOCATE(gradient_i(src_size), gradient_j(src_size), & 00415 gradient_ij(src_size)) 00416 00417 DO ib=1, nbr_comp 00418 CALL gradient_bicubic(nlon_src,nlat_src,src_array(:,ib), & 00419 sou_mask, src_lat, src_lon, id_sper, cd_sper, & 00420 gradient_i, gradient_j, gradient_ij) 00421 00422 DO n = 1, num_links 00423 IF (src_addr(n) .NE. 0) THEN 00424 dst_array(dst_addr(n),ib) = dst_array(dst_addr(n),ib) & 00425 + weights(1,n) * src_array(src_addr(n),ib) & 00426 + weights(2,n) * gradient_i(src_addr(n)) & 00427 + weights(3,n) * gradient_j(src_addr(n)) & 00428 + weights(4,n) * gradient_ij(src_addr(n)) 00429 weightot(dst_addr(n),ib) = weightot(dst_addr(n),ib) & 00430 + weights(1,n) + weights(2,n) + weights(3,n) + weights(4,n) 00431 ENDIF 00432 ENDDO 00433 END DO 00434 00435 DEALLOCATE(gradient_i, gradient_j, gradient_ij) 00436 CASE ('D') !reduced 00437 DO ib=1, nbr_comp 00438 DO n = 1, num_links 00439 IF (src_addr(n) .NE. 0) THEN 00440 dst_array(dst_addr(n),ib) = dst_array(dst_addr(n),ib) & 00441 + weights(1,n) * src_array(src_addr(n),ib) 00442 weightot(dst_addr(n),ib) = weightot(dst_addr(n),ib) & 00443 + weights(1,n) 00444 ENDIF 00445 ENDDO 00446 ENDDO 00447 END SELECT 00448 CASE ('DISTWGT') ! distance weighted average 00449 DO ib=1, nbr_comp 00450 DO n = 1, num_links 00451 IF (src_addr(n) .NE. 0) THEN 00452 dst_array(dst_addr(n),ib) = dst_array(dst_addr(n),ib) & 00453 + weights(1,n) * src_array(src_addr(n),ib) 00454 weightot(dst_addr(n),ib) = weightot(dst_addr(n),ib) & 00455 + weights(1,n) 00456 ENDIF 00457 ENDDO 00458 ENDDO 00459 CASE ('GAUSWGT') ! distance gaussian weighted average 00460 DO ib=1, nbr_comp 00461 DO n = 1, num_links 00462 IF (src_addr(n) .NE. 0) THEN 00463 dst_array(dst_addr(n),ib) = dst_array(dst_addr(n),ib) & 00464 + weights(1,n) * src_array(src_addr(n),ib) 00465 weightot(dst_addr(n),ib) = weightot(dst_addr(n),ib) & 00466 + weights(1,n) 00467 ENDIF 00468 ENDDO 00469 ENDDO 00470 END SELECT ! remapping method 00471 00472 IF (ll_weightot) THEN 00473 DO ib=1, nbr_comp 00474 DO n = 1, dst_size 00475 IF (weightot(n,ib) .LT. EPSILON(1.)) dst_array(n,ib) = 1.0E+20 00476 ENDDO 00477 ENDDO 00478 ENDIF 00479 00480 DEALLOCATE (src_addr, dst_addr, weights) 00481 00482 ! 00483 IF (nlogprt .GE. 2) THEN 00484 WRITE (UNIT = nulou,FMT = *)' ' 00485 WRITE (UNIT = nulou,FMT = *)'Leaving routine remap_vector_comps' 00486 CALL FLUSH(nulou) 00487 ENDIF 00488 ! 00489 END SUBROUTINE remap_vector_comps 00490 00491 ! 00492 !** ----------------------------------------------------------------------- 00493 !** ----------------------------------------------------------------------- 00494 ! 00495 00496 SUBROUTINE check_points_at_poles (arrayA, arrayB, latA, latB, grd_nameA, & 00497 grd_nameB, dst_size) 00498 ! 00499 !** Description : 00500 ! ------------- 00501 ! Test if there are points on the latitudes 90deg north and 90deg south. 00502 ! If there are points on these latitudes, the averege for these points 00503 ! in arrayA and arrayB are calculated and distributed. arrayA and arrayB 00504 ! are modified if necessary. 00505 ! 00506 IMPLICIT NONE 00507 ! 00508 !** Input variables 00509 ! ---------------- 00510 INTEGER(KIND=int_kind), INTENT(in) :: dst_size 00511 REAL (KIND=real_kind), DIMENSION(dst_size), INTENT(in):: latA, latB 00512 CHARACTER(LEN=8), INTENT(in) :: grd_nameA, grd_nameB 00513 ! 00514 !** OutInput variables 00515 ! ------------------ 00516 REAL (KIND=real_kind), DIMENSION(dst_size), INTENT(inout):: 00517 arrayA, ! VECTOR_I field on target grid dstA 00518 arrayB ! VECTOR_J field on target grid dstB 00519 ! 00520 !** Local variables 00521 ! ---------------- 00522 INTEGER(KIND=int_kind) :: nbr_N, nbr_S, n 00523 REAL(kind=real_kind) :: ave_N, ave_S, latpolN, latpolS 00524 ! 00525 !** ---------------------------------------------------------------------------- 00526 ! 00527 IF (nlogprt .GE. 2) THEN 00528 WRITE (UNIT = nulou,FMT = *)' ' 00529 WRITE (UNIT = nulou,FMT = *)'Entering routine check_points_at_poles' 00530 CALL FLUSH(nulou) 00531 ENDIF 00532 ! 00533 latpolN = pi*half 00534 latpolS = -pi*half 00535 nbr_N = 0 00536 nbr_S = 0 00537 ave_N = 0.0 00538 ave_S = 0.0 00539 00540 DO n = 1, dst_size 00541 IF (latA(n) == latpolN) THEN 00542 nbr_N = nbr_N + 1 00543 ave_N= ave_N + arrayA(n) 00544 ELSE IF (latA(n) == latpolS) THEN 00545 nbr_S = nbr_S + 1 00546 ave_S= ave_S + arrayA(n) 00547 END IF 00548 END DO 00549 00550 IF (nbr_N .NE. 0) THEN 00551 ave_N = ave_N/nbr_N 00552 WHERE (latA == latpolN) 00553 arrayA = ave_N 00554 END WHERE 00555 END IF 00556 00557 IF (nbr_S .NE. 0) THEN 00558 ave_S = ave_S/nbr_S 00559 WHERE (latA .EQ. latpolS) 00560 arrayA = ave_S 00561 END WHERE 00562 END IF 00563 ! 00564 IF (nlogprt .GE. 2) THEN 00565 WRITE (UNIT = nulou,FMT = *) ' ' 00566 WRITE (UNIT = nulou,FMT = *) 'For target grid : ', grd_nameA 00567 WRITE (UNIT = nulou,FMT = *) nbr_N,' points at the north pole ' 00568 WRITE (UNIT = nulou,FMT = *) nbr_S,' points at the south pole ' 00569 WRITE (UNIT = nulou,FMT = *) ' ' 00570 IF (nbr_N/=0 .OR. nbr_S/=0) THEN 00571 WRITE (UNIT = nulou,FMT = *) 'Average of field component I at north pole : ', ave_N 00572 WRITE (UNIT = nulou,FMT = *) 'Average of field component I at south pole : ', ave_S 00573 ENDIF 00574 CALL FLUSH (nulou) 00575 ENDIF 00576 ! 00577 !** The same calculation for the second target grid 00578 ! 00579 nbr_N = 0 00580 nbr_S = 0 00581 ave_N = 0.0 00582 ave_S = 0.0 00583 00584 DO n = 1, dst_size 00585 IF (latB(n) == latpolN) THEN 00586 nbr_N = nbr_N + 1 00587 ave_N= ave_N + arrayB(n) 00588 ELSE IF (latB(n) == latpolS) THEN 00589 nbr_S = nbr_S + 1 00590 ave_S= ave_S + arrayB(n) 00591 END IF 00592 END DO 00593 00594 IF (nbr_N .NE. 0) THEN 00595 ave_N = ave_N/nbr_N 00596 WHERE (latB == latpolN) 00597 arrayB = ave_N 00598 END WHERE 00599 END IF 00600 00601 IF (nbr_S .NE. 0) THEN 00602 ave_S = ave_S/nbr_S 00603 WHERE (latB .EQ. latpolS) 00604 arrayB = ave_S 00605 END WHERE 00606 END IF 00607 ! 00608 IF (nlogprt .GE. 2) THEN 00609 IF (grd_nameA .NE. grd_nameB) THEN 00610 WRITE (UNIT = nulou,FMT = *) ' ' 00611 WRITE (UNIT = nulou,FMT = *) 'For target grid : ', grd_nameB 00612 WRITE (UNIT = nulou,FMT = *) nbr_N,' points at the north pole ' 00613 WRITE (UNIT = nulou,FMT = *) nbr_S,' points at the south pole ' 00614 WRITE (UNIT = nulou,FMT = *) ' ' 00615 ENDIF 00616 IF (nbr_N/=0 .OR. nbr_S/=0) THEN 00617 WRITE (UNIT = nulou,FMT = *) 'Average of field component J at north pole : ', ave_N 00618 WRITE (UNIT = nulou,FMT = *) 'Average of field component J at south pole : ', ave_S 00619 ENDIF 00620 ! 00621 WRITE (UNIT = nulou,FMT = *)' ' 00622 WRITE (UNIT = nulou,FMT = *)'Leaving routine check_points_at_poles' 00623 CALL FLUSH(nulou) 00624 ENDIF 00625 ! 00626 END SUBROUTINE check_points_at_poles 00627 00628 ! 00629 !** ----------------------------------------------------------------------- 00630 !** ----------------------------------------------------------------------- 00631 ! 00632 00633 SUBROUTINE write_src_array_spheric(array, array_size, nlon, nlat, & 00634 grd_name_src, grd_name_dst) 00635 ! 00636 !** Description 00637 ! ----------- 00638 ! This routine will write the 2 spherical components given in 'array' in 00639 ! netcdf file vector_debug.nc. The file will first be created. 00640 ! 00641 IMPLICIT NONE 00642 ! 00643 !** Input variables 00644 ! 00645 INTEGER (KIND=int_kind), INTENT(in) :: array_size, nlon, nlat 00646 REAL (KIND=real_kind), DIMENSION(array_size,2), INTENT(in) :: array 00647 CHARACTER(LEN=8), INTENT(in) :: grd_name_src, grd_name_dst 00648 ! 00649 !** Local variables 00650 ! 00651 INTEGER(KIND=int_kind) :: dim_id(2), nc_fileid, dimI_id, dimJ_id, stat, 00652 varI_id, varJ_id, icount, kcount 00653 INTEGER(KIND=int_kind) :: ilenstr 00654 CHARACTER(char_len) :: text, filename 00655 ! 00656 !** -------------------------------------------------------------------------- 00657 ! 00658 IF (nlogprt .GE. 2) THEN 00659 WRITE(nulou,*) ' ' 00660 WRITE(nulou,*) ' Entering routine write_src_array_spheric' 00661 call flush(nulou) 00662 ENDIF 00663 ! 00664 ! 00665 !** Create the file vector_debug.nc and define the dimensions 00666 ! 00667 00668 icount = ilenstr(grd_name_src,jpeight) 00669 kcount = ilenstr(grd_name_dst, jpeight) 00670 filename='vector_debug_'//grd_name_src(1:icount)//'_to_'//grd_name_dst(1:kcount)//'.nc' 00671 CALL hdlerr(NF_CREATE(filename, 0, nc_fileid), 'write_src_array_spheric') 00672 ! 00673 IF (nlogprt .GE. 2) THEN 00674 WRITE (UNIT = nulou,FMT = *) ' ' 00675 WRITE (UNIT = nulou,FMT = *) filename,' is created containg fields in spheric and' 00676 WRITE (UNIT = nulou,FMT = *) ' eventually cartesian referentials' 00677 WRITE (UNIT = nulou,FMT = *) ' ' 00678 CALL FLUSH(nulou) 00679 ENDIF 00680 ! 00681 CALL hdlerr(NF_DEF_DIM & 00682 (nc_fileid,'src_grid_dim_I',nlon,dimI_id), 'write_src_array_spheric') 00683 00684 CALL hdlerr(NF_DEF_DIM & 00685 (nc_fileid,'src_grid_dim_J',nlat,dimJ_id), 'write_src_array_spheric') 00686 00687 dim_id(1) = dimI_id 00688 dim_id(2) = dimJ_id 00689 ! 00690 !** Define variable for component I 00691 ! 00692 CALL hdlerr(NF_DEF_VAR & 00693 (nc_fileid,'Comp_I_src_spheric',NF_DOUBLE,2,dim_id,varI_id), & 00694 'write_src_array_spheric') 00695 00696 text = 'Component I in spherical referential on source grid '//grd_name_src 00697 00698 CALL hdlerr(NF_PUT_ATT_TEXT(nc_fileid,varI_id, 'array',& 00699 LEN(text),text),'write_src_array_spheric') 00700 ! 00701 !** Define variable for component J 00702 ! 00703 CALL hdlerr(NF_DEF_VAR & 00704 (nc_fileid,'Comp_J_src_spheric',NF_DOUBLE,2,dim_id,varJ_id), & 00705 'write_src_array_spheric') 00706 00707 text = 'Component J in spherical referential on source grid '//grd_name_src 00708 00709 CALL hdlerr(NF_PUT_ATT_TEXT(nc_fileid,varJ_id, 'array',& 00710 LEN(text),text),'write_src_array_spheric') 00711 00712 CALL hdlerr(NF_ENDDEF(nc_fileid), 'write_src_array_spheric') 00713 00714 ! 00715 !** Put the component I and J into file 00716 ! 00717 CALL hdlerr(NF_PUT_VAR_DOUBLE & 00718 (nc_fileid,varI_id,array(:,1)),'write_src_array_spheric') 00719 00720 CALL hdlerr(NF_PUT_VAR_DOUBLE & 00721 (nc_fileid,varJ_id,array(:,2)),'write_src_array_spheric') 00722 00723 CALL hdlerr(NF_CLOSE(nc_fileid),'write_src_array_spheric') 00724 ! 00725 IF (nlogprt .GE. 2) THEN 00726 WRITE(nulou,*) ' ' 00727 WRITE(nulou,*) ' Entering routine write_src_array_spheric' 00728 call flush(nulou) 00729 ENDIF 00730 ! 00731 END SUBROUTINE write_src_array_spheric 00732 00733 ! 00734 !** ----------------------------------------------------------------------- 00735 !** ----------------------------------------------------------------------- 00736 ! 00737 00738 SUBROUTINE write_dst_array_spheric(array, array_size, nlon, nlat, & 00739 grd_name_src, grd_name_dst) 00740 ! 00741 !** Description 00742 ! ----------- 00743 ! This routine will write the 3 spherical components given in 'array' in 00744 ! netcdf file vector_debug_XXX_to_YYY.nc. The file has to existe(it is created by 00745 ! call to write_src_array_spheric). 00746 ! 00747 IMPLICIT NONE 00748 ! 00749 !** Input variables 00750 ! --------------- 00751 INTEGER (KIND=int_kind), INTENT(in) :: array_size, nlon, nlat 00752 REAL (KIND=real_kind), DIMENSION(array_size,3), INTENT(in) :: array 00753 CHARACTER(LEN=8), INTENT(in) :: grd_name_src, grd_name_dst 00754 ! 00755 !** Local variables 00756 ! --------------- 00757 INTEGER(KIND=int_kind) :: dim_id(2), nc_fileid, dimI_id, dimJ_id, 00758 stat, varI_id, varJ_id, varK_id, icount, kcount 00759 INTEGER(KIND=int_kind) :: ilenstr 00760 CHARACTER(char_len) :: text, filename 00761 ! 00762 !** ----------------------------------------------------------------------- 00763 ! 00764 IF (nlogprt .GE. 2) THEN 00765 WRITE(nulou,*) ' ' 00766 WRITE(nulou,*) ' Entering routine write_dst_array_spheric' 00767 call flush(nulou) 00768 ENDIF 00769 ! 00770 !** Open the file vector_debug_XXX_to_YYY.nc and define the dimensions 00771 ! 00772 icount = ilenstr(grd_name_src,jpeight) 00773 kcount = ilenstr(grd_name_dst, jpeight) 00774 filename='vector_debug_'//grd_name_src(1:icount)//'_to_'//grd_name_dst(1:kcount)//'.nc' 00775 CALL hdlerr(NF_OPEN(filename, NF_WRITE, nc_fileid), 'write_dst_array_spheric') 00776 00777 CALL hdlerr(NF_REDEF(nc_fileid), 'write_dst_array_spheric') 00778 00779 stat=NF_INQ_DIMID(nc_fileid,'dst_grid_dim_I',dimI_id) 00780 IF (stat .NE. NF_NOERR) CALL hdlerr(NF_DEF_DIM & 00781 (nc_fileid,'dst_grid_dim_I',nlon,dimI_id), 'write_dst_array_spheric') 00782 00783 stat=NF_INQ_DIMID(nc_fileid,'dst_grid_dim_J',dimJ_id) 00784 IF (stat .NE. NF_NOERR) CALL hdlerr(NF_DEF_DIM & 00785 (nc_fileid,'dst_grid_dim_J',nlat,dimJ_id), 'write_dst_array_spheric') 00786 00787 dim_id(1) = dimI_id 00788 dim_id(2) = dimJ_id 00789 ! 00790 !** Define variable for component I 00791 ! 00792 CALL hdlerr(NF_DEF_VAR & 00793 (nc_fileid,'Comp_I_dst_spheric',NF_DOUBLE,2,dim_id,varI_id), & 00794 'write_dst_array_spheric') 00795 00796 text = 'Component I in spherical referential on target grid '//grd_name_dst 00797 00798 CALL hdlerr(NF_PUT_ATT_TEXT(nc_fileid,varI_id, 'array',& 00799 LEN(text),text),'write_dst_array_spheric') 00800 ! 00801 !* Defining variable for component J 00802 ! 00803 CALL hdlerr(NF_DEF_VAR & 00804 (nc_fileid,'Comp_J_dst_spheric',NF_DOUBLE,2,dim_id,varJ_id), & 00805 'write_dst_array_spheric') 00806 00807 text = 'Component J in spherical referential on target grid '//grd_name_dst 00808 CALL hdlerr(NF_PUT_ATT_TEXT(nc_fileid,varJ_id, 'array',& 00809 LEN(text),text),'write_dst_array_spheric') 00810 ! 00811 !* Defining variable for component K 00812 ! 00813 CALL hdlerr(NF_DEF_VAR & 00814 (nc_fileid,'Comp_K_dst_spheric',NF_DOUBLE,2,dim_id,varK_id), & 00815 'write_dst_array_spheric') 00816 00817 text = 'Component K in spherical referential on target grid '//grd_name_dst 00818 CALL hdlerr(NF_PUT_ATT_TEXT(nc_fileid,varK_id, 'array',& 00819 LEN(text),text),'write_dst_array_spheric') 00820 00821 CALL hdlerr(NF_ENDDEF(nc_fileid), 'write_src_array_spheric') 00822 ! 00823 !* Putting the fields into the file 00824 ! 00825 CALL hdlerr(NF_PUT_VAR_DOUBLE & 00826 (nc_fileid,varI_id,array(:,1)),'write_dst_array_spheric') 00827 00828 CALL hdlerr(NF_PUT_VAR_DOUBLE & 00829 (nc_fileid,varJ_id,array(:,2)),'write_dst_array_spheric') 00830 00831 CALL hdlerr(NF_PUT_VAR_DOUBLE & 00832 (nc_fileid,varK_id,array(:,3)),'write_dst_array_spheric') 00833 00834 CALL hdlerr(NF_CLOSE(nc_fileid),'write_src_array_spheric') 00835 ! 00836 IF (nlogprt .GE. 2) THEN 00837 WRITE(nulou,*) ' ' 00838 WRITE(nulou,*) ' Leaving routine write_dst_array_spheric' 00839 call flush(nulou) 00840 ENDIF 00841 ! 00842 END SUBROUTINE write_dst_array_spheric 00843 00844 ! 00845 !** ----------------------------------------------------------------------- 00846 !** ----------------------------------------------------------------------- 00847 ! 00848 00849 SUBROUTINE write_cartesian_components(src_array, dst_array, src_size, dst_size, & 00850 nlon_src, nlat_src, nlon_dst, nlat_dst, grd_name_src, grd_name_dst) 00851 ! 00852 !** Description 00853 ! ----------- 00854 ! This routine will write the 3 cartesian components on soure and target grid in 00855 ! netcdf file vector_debug_XXX_to_YYY.nc. The file has to existe(it is created by 00856 ! call to write_src_array_spheric). 00857 ! 00858 IMPLICIT NONE 00859 ! 00860 !** Input variables 00861 ! --------------- 00862 INTEGER (KIND=int_kind), INTENT(in) :: src_size, dst_size, 00863 nlon_src, nlat_src, nlon_dst, nlat_dst 00864 REAL (KIND=real_kind), DIMENSION(src_size,3), INTENT(in) :: src_array 00865 REAL (KIND=real_kind), DIMENSION(dst_size,3), INTENT(in) :: dst_array 00866 CHARACTER(LEN=8), INTENT(in) :: grd_name_src, grd_name_dst 00867 ! 00868 !** Local variables 00869 ! --------------- 00870 ! 00871 INTEGER(KIND=int_kind) :: dimsrc_id(2), dimdst_id(2), nc_fileid, 00872 dimI_id, dimJ_id, stat, 00873 varXsrc_id, varXdst_id, varYsrc_id, varYdst_id, varZsrc_id, varZdst_id, 00874 icount, kcount 00875 INTEGER(KIND=int_kind) :: ilenstr 00876 CHARACTER(char_len) :: text, filename 00877 ! 00878 !** ---------------------------------------------------------------------------- 00879 ! 00880 IF (nlogprt .GE. 2) THEN 00881 WRITE(nulou,*) ' ' 00882 WRITE(nulou,*) ' Entering routine write_cartesian_components' 00883 call flush(nulou) 00884 ENDIF 00885 ! 00886 !** Open the file vector_debug_XXX_to_YYY.nc and define the dimensions 00887 ! 00888 00889 icount = ilenstr(grd_name_src,jpeight) 00890 kcount = ilenstr(grd_name_dst, jpeight) 00891 filename='vector_debug_'//grd_name_src(1:icount)//'_to_'//grd_name_dst(1:kcount)//'.nc' 00892 CALL hdlerr(NF_OPEN(filename, NF_WRITE, nc_fileid), 'write_cartesian_components') 00893 00894 CALL hdlerr(NF_REDEF(nc_fileid), 'write_cartesian_components') 00895 ! 00896 !** Define dimension for src grid 00897 ! 00898 stat=NF_INQ_DIMID(nc_fileid,'src_grid_dim_I',dimI_id) 00899 IF (stat .NE. NF_NOERR) CALL hdlerr(NF_DEF_DIM & 00900 (nc_fileid,'src_grid_dim_I',nlon_src,dimI_id), 'write_cartesian_components') 00901 00902 stat=NF_INQ_DIMID(nc_fileid,'src_grid_dim_J',dimJ_id) 00903 IF (stat .NE. NF_NOERR) CALL hdlerr(NF_DEF_DIM & 00904 (nc_fileid,'src_grid_dim_J',nlat_src,dimJ_id), 'write_cartesian_components') 00905 00906 dimsrc_id(1) = dimI_id 00907 dimsrc_id(2) = dimJ_id 00908 ! 00909 !** Define dimension for dst grid 00910 ! 00911 stat=NF_INQ_DIMID(nc_fileid,'dst_grid_dim_I',dimI_id) 00912 IF (stat .NE. NF_NOERR) CALL hdlerr(NF_DEF_DIM & 00913 (nc_fileid,'dst_grid_dim_I',nlon_dst,dimI_id), 'write_cartesian_components') 00914 00915 stat=NF_INQ_DIMID(nc_fileid,'dst_grid_dim_J',dimJ_id) 00916 IF (stat .NE. NF_NOERR) CALL hdlerr(NF_DEF_DIM & 00917 (nc_fileid,'dst_grid_dim_J',nlat_dst,dimJ_id), 'write_cartesian_components') 00918 00919 dimdst_id(1) = dimI_id 00920 dimdst_id(2) = dimJ_id 00921 00922 ! 00923 !* Define src and dst variables for component X 00924 ! 00925 CALL hdlerr(NF_DEF_VAR & 00926 (nc_fileid,'Comp_X_cart_src',NF_DOUBLE,2,dimsrc_id, varXsrc_id), & 00927 'write_cartesian_components') 00928 text = 'Component X in cartesian referential on source grid '//grd_name_src 00929 CALL hdlerr(NF_PUT_ATT_TEXT(nc_fileid,varXsrc_id, 'array',& 00930 LEN(text),text),'write_cartesian_components') 00931 00932 CALL hdlerr(NF_DEF_VAR & 00933 (nc_fileid,'Comp_X_cart_dst',NF_DOUBLE,2,dimdst_id, varXdst_id), & 00934 'write_cartesian_components') 00935 text = 'Component X in cartesian referential on target grid '//grd_name_dst 00936 CALL hdlerr(NF_PUT_ATT_TEXT(nc_fileid,varXdst_id, 'array',& 00937 LEN(text),text),'write_cartesian_components') 00938 00939 00940 ! 00941 !* Define src and dst variables for component Y 00942 ! 00943 CALL hdlerr(NF_DEF_VAR & 00944 (nc_fileid,'Comp_Y_cart_src',NF_DOUBLE,2,dimsrc_id, varYsrc_id), & 00945 'write_cartesian_components') 00946 text = 'Component Y in cartesian referential on source grid '//grd_name_src 00947 CALL hdlerr(NF_PUT_ATT_TEXT(nc_fileid,varYsrc_id, 'array',& 00948 LEN(text),text),'write_cartesian_components') 00949 00950 CALL hdlerr(NF_DEF_VAR & 00951 (nc_fileid,'Comp_Y_cart_dst',NF_DOUBLE,2,dimdst_id, varYdst_id), & 00952 'write_cartesian_components') 00953 text = 'Component Y in cartesian referential on target grid '//grd_name_dst 00954 CALL hdlerr(NF_PUT_ATT_TEXT(nc_fileid,varYdst_id, 'array',& 00955 LEN(text),text),'write_cartesian_components') 00956 00957 00958 ! 00959 !* Define src and dst variables for component Z 00960 ! 00961 CALL hdlerr(NF_DEF_VAR & 00962 (nc_fileid,'Comp_Z_cart_src',NF_DOUBLE,2,dimsrc_id, varZsrc_id), & 00963 'write_cartesian_components') 00964 text = 'Component Z in cartesian referential on source grid '//grd_name_src 00965 CALL hdlerr(NF_PUT_ATT_TEXT(nc_fileid,varZsrc_id, 'array',& 00966 LEN(text),text),'write_cartesian_components') 00967 00968 CALL hdlerr(NF_DEF_VAR & 00969 (nc_fileid,'Comp_Z_cart_dst',NF_DOUBLE,2,dimdst_id, varZdst_id), & 00970 'write_cartesian_components') 00971 text = 'Component Z in cartesian referential on target grid '//grd_name_dst 00972 CALL hdlerr(NF_PUT_ATT_TEXT(nc_fileid,varZdst_id, 'array',& 00973 LEN(text),text),'write_cartesian_components') 00974 00975 00976 CALL hdlerr(NF_ENDDEF(nc_fileid), 'write_src_array_spheric') 00977 00978 ! 00979 !* Putting the fields into the file 00980 ! 00981 CALL hdlerr(NF_PUT_VAR_DOUBLE & 00982 (nc_fileid,varXsrc_id,src_array(:,1)),'write_cartesian_components') 00983 CALL hdlerr(NF_PUT_VAR_DOUBLE & 00984 (nc_fileid,varXdst_id,dst_array(:,1)),'write_cartesian_components') 00985 00986 CALL hdlerr(NF_PUT_VAR_DOUBLE & 00987 (nc_fileid,varYsrc_id,src_array(:,2)),'write_cartesian_components') 00988 CALL hdlerr(NF_PUT_VAR_DOUBLE & 00989 (nc_fileid,varYdst_id,dst_array(:,2)),'write_cartesian_components') 00990 00991 CALL hdlerr(NF_PUT_VAR_DOUBLE & 00992 (nc_fileid,varZsrc_id,src_array(:,3)),'write_cartesian_components') 00993 CALL hdlerr(NF_PUT_VAR_DOUBLE & 00994 (nc_fileid,varZdst_id,dst_array(:,3)),'write_cartesian_components') 00995 00996 00997 CALL hdlerr(NF_CLOSE(nc_fileid),'write_src_array_spheric') 00998 ! 00999 IF (nlogprt .GE. 2) THEN 01000 WRITE(nulou,*) ' ' 01001 WRITE(nulou,*) ' Leaving routine write_cartesian_components' 01002 call flush(nulou) 01003 ENDIF 01004 ! 01005 END SUBROUTINE write_cartesian_components 01006 01007 ! 01008 !** ----------------------------------------------------------------------- 01009 !** ----------------------------------------------------------------------- 01010 ! 01011 01012 END MODULE vector
1.7.4 
