- Timestamp:
- 2020-05-05T08:18:05+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/SBC/fldread.F90
r12738 r12866 71 71 CHARACTER(len = 8) :: cltype ! type of data file 'daily', 'monthly' or yearly' 72 72 INTEGER :: num ! iom id of the jpfld files to be read 73 INTEGER , DIMENSION(2) :: nrec_b ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year) 74 INTEGER , DIMENSION(2) :: nrec_a ! after record (1: index, 2: second since Jan. 1st 00h of nit000 year) 73 INTEGER , DIMENSION(2,2) :: nrec ! before/after record (1: index, 2: second since Jan. 1st 00h of yr nit000) 74 INTEGER :: nbb ! index of before values 75 INTEGER :: naa ! index of after values 75 76 INTEGER , ALLOCATABLE, DIMENSION(: ) :: nrecsec ! 76 77 REAL(wp), ALLOCATABLE, DIMENSION(:,:,: ) :: fnow ! input fields interpolated to now time step … … 156 157 INTEGER :: jf ! dummy indices 157 158 INTEGER :: isecsbc ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step 159 INTEGER :: ibb, iaa ! shorter name for sd(jf)%nbb and sd(jf)%naa 158 160 LOGICAL :: ll_firstcall ! true if this is the first call to fld_read for this set of fields 159 161 REAL(wp) :: zt_offset ! local time offset variable … … 203 205 IF( TRIM(sd(jf)%clrootname) == 'NOT USED' ) CYCLE 204 206 ! 207 ibb = sd(jf)%nbb ; iaa = sd(jf)%naa 208 ! 205 209 IF( sd(jf)%ln_tint ) THEN ! temporal interpolation 206 210 IF(lwp .AND. kt - nit000 <= 100 ) THEN … … 208 212 & "', records b/a: ', i6.4, '/', i6.4, ' (days ', f9.4,'/', f9.4, ')')" 209 213 WRITE(numout, clfmt) TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday, & 210 & sd(jf)%nrec _b(1), sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday214 & sd(jf)%nrec(1,ibb), sd(jf)%nrec(1,iaa), REAL(sd(jf)%nrec(2,ibb),wp)/rday, REAL(sd(jf)%nrec(2,iaa),wp)/rday 211 215 WRITE(numout, *) ' zt_offset is : ',zt_offset 212 216 ENDIF 213 217 ! temporal interpolation weights 214 ztinta = REAL( isecsbc - sd(jf)%nrec _b(2), wp ) / REAL( sd(jf)%nrec_a(2) - sd(jf)%nrec_b(2), wp )218 ztinta = REAL( isecsbc - sd(jf)%nrec(2,ibb), wp ) / REAL( sd(jf)%nrec(2,iaa) - sd(jf)%nrec(2,ibb), wp ) 215 219 ztintb = 1. - ztinta 216 sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:, 1) + ztinta * sd(jf)%fdta(:,:,:,2)220 sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:,ibb) + ztinta * sd(jf)%fdta(:,:,:,iaa) 217 221 ELSE ! nothing to do... 218 222 IF(lwp .AND. kt - nit000 <= 100 ) THEN … … 220 224 & "', record: ', i6.4, ' (days ', f9.4, ' <-> ', f9.4, ')')" 221 225 WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday, & 222 & sd(jf)%nrec _a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday226 & sd(jf)%nrec(1,iaa), REAL(sd(jf)%nrec(2,ibb),wp)/rday, REAL(sd(jf)%nrec(2,iaa),wp)/rday 223 227 ENDIF 224 228 ENDIF … … 250 254 ! 251 255 CALL fld_clopn( sdjf ) 252 sdjf%nrec _a(:) = (/ 1, nflag /) ! default definition to force flp_update to read the file.256 sdjf%nrec(:,sdjf%naa) = (/ 1, nflag /) ! default definition to force flp_update to read the file. 253 257 ! 254 258 END SUBROUTINE fld_init … … 261 265 !! ** Purpose : Compute 262 266 !! if sdjf%ln_tint = .TRUE. 263 !! nrec _a: record number and its time (nrec_b is obtained from nrec_awhen swapping)267 !! nrec(:,iaa): record number and its time (nrec(:,ibb) is obtained from nrec(:,iaa) when swapping) 264 268 !! if sdjf%ln_tint = .FALSE. 265 !! nrec _a(1): record number266 !! nrec _b(2) and nrec_a(2): time of the beginning and end of the record269 !! nrec(1,iaa): record number 270 !! nrec(2,ibb) and nrec(2,iaa): time of the beginning and end of the record 267 271 !!---------------------------------------------------------------------- 268 272 INTEGER , INTENT(in ) :: ksecsbc ! … … 270 274 INTEGER , OPTIONAL, INTENT(in ) :: Kmm ! ocean time level index 271 275 ! 272 INTEGER :: ja ! end of this record (in seconds) 273 !!---------------------------------------------------------------------- 274 ! 275 IF( ksecsbc > sdjf%nrec_a(2) ) THEN ! --> we need to update after data 276 INTEGER :: ja ! end of this record (in seconds) 277 INTEGER :: ibb, iaa ! shorter name for sdjf%nbb and sdjf%naa 278 !!---------------------------------------------------------------------- 279 ibb = sdjf%nbb ; iaa = sdjf%naa 280 ! 281 IF( ksecsbc > sdjf%nrec(2,iaa) ) THEN ! --> we need to update after data 276 282 277 ! find where is the new after record... (it is not necessary sdjf%nrec _a(1)+1 )278 ja = sdjf%nrec _a(1)283 ! find where is the new after record... (it is not necessary sdjf%nrec(1,iaa)+1 ) 284 ja = sdjf%nrec(1,iaa) 279 285 DO WHILE ( ksecsbc >= sdjf%nrecsec(ja) .AND. ja < sdjf%nreclast ) ! Warning: make sure ja <= sdjf%nreclast in this test 280 286 ja = ja + 1 … … 283 289 284 290 ! if ln_tint and if the new after is not ja+1, we need also to update after data before the swap 285 ! so, after the swap, sdjf%nrec _b(2) will still be the closest value located just before ksecsbc286 IF( sdjf%ln_tint .AND. ( ja > sdjf%nrec _a(1) + 1 .OR. sdjf%nrec_a(2) == nflag ) ) THEN287 sdjf%nrec _a(:) = (/ ja-1, sdjf%nrecsec(ja-1) /) ! update nrec_awith before information288 CALL fld_get( sdjf, Kmm ) ! read after data that will be used as before data291 ! so, after the swap, sdjf%nrec(2,ibb) will still be the closest value located just before ksecsbc 292 IF( sdjf%ln_tint .AND. ( ja > sdjf%nrec(1,iaa) + 1 .OR. sdjf%nrec(2,iaa) == nflag ) ) THEN 293 sdjf%nrec(:,iaa) = (/ ja-1, sdjf%nrecsec(ja-1) /) ! update nrec(:,iaa) with before information 294 CALL fld_get( sdjf, Kmm ) ! read after data that will be used as before data 289 295 ENDIF 290 296 … … 309 315 ! if ln_tint and if after is not the first record, we must (potentially again) update after data before the swap 310 316 IF( sdjf%ln_tint .AND. ja > 1 ) THEN 311 IF( sdjf%nrecsec(0) /= nflag ) THEN ! no trick used: after file is not the current file312 sdjf%nrec _a(:) = (/ ja-1, sdjf%nrecsec(ja-1) /) ! update nrec_awith before information313 CALL fld_get( sdjf, Kmm ) ! read after data that will be used as before data317 IF( sdjf%nrecsec(0) /= nflag ) THEN ! no trick used: after file is not the current file 318 sdjf%nrec(:,iaa) = (/ ja-1, sdjf%nrecsec(ja-1) /) ! update nrec(:,iaa) with before information 319 CALL fld_get( sdjf, Kmm ) ! read after data that will be used as before data 314 320 ENDIF 315 321 ENDIF … … 317 323 ENDIF 318 324 319 IF( sdjf%ln_tint ) THEN 320 ! Swap data 321 sdjf%nrec_b(:) = sdjf%nrec_a(:) ! swap before record informations 322 sdjf%rotn(1) = sdjf%rotn(2) ! swap before rotate informations 323 sdjf%fdta(:,:,:,1) = sdjf%fdta(:,:,:,2) ! swap before record field 324 ELSE 325 sdjf%nrec_b(:) = (/ ja-1, sdjf%nrecsec(ja-1) /) ! only for print 325 IF( sdjf%ln_tint ) THEN ! Swap data 326 sdjf%nbb = sdjf%naa ! swap indices 327 sdjf%naa = 3 - sdjf%naa ! = 2(1) if naa == 1(2) 328 ELSE ! No swap 329 sdjf%nrec(:,ibb) = (/ ja-1, sdjf%nrecsec(ja-1) /) ! only for print 326 330 ENDIF 327 331 328 332 ! read new after data 329 sdjf%nrec _a(:) = (/ ja, sdjf%nrecsec(ja) /) ! update nrec_aas it is used by fld_get330 CALL fld_get( sdjf, Kmm ) ! read after data (with nrec_ainformations)333 sdjf%nrec(:,sdjf%naa) = (/ ja, sdjf%nrecsec(ja) /) ! update nrec(:,naa) as it is used by fld_get 334 CALL fld_get( sdjf, Kmm ) ! read after data (with nrec(:,naa) informations) 331 335 332 336 ENDIF … … 345 349 ! 346 350 INTEGER :: ipk ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) 351 INTEGER :: iaa ! shorter name for sdjf%naa 347 352 INTEGER :: iw ! index into wgts array 348 353 INTEGER :: ipdom ! index of the domain … … 353 358 ! 354 359 ipk = SIZE( sdjf%fnow, 3 ) 360 iaa = sdjf%naa 355 361 ! 356 362 IF( ASSOCIATED(sdjf%imap) ) THEN 357 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:, 2), sdjf%nrec_a(1), &363 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,iaa), sdjf%nrec(1,iaa), & 358 364 & sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel, sdjf%lzint, Kmm ) 359 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1), &365 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec(1,iaa), & 360 366 & sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel, sdjf%lzint, Kmm ) 361 367 ENDIF 362 368 ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 363 369 CALL wgt_list( sdjf, iw ) 364 IF( sdjf%ln_tint ) THEN ; CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fdta(:,:,:,2), & 365 & sdjf%nrec_a(1), sdjf%lsmname ) 366 ELSE ; CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fnow(:,:,: ), & 367 & sdjf%nrec_a(1), sdjf%lsmname ) 370 IF( sdjf%ln_tint ) THEN ; CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fdta(:,:,:,iaa), & 371 & sdjf%nrec(1,iaa), sdjf%lsmname ) 372 CALL lbc_lnk( 'fldread', sdjf%fdta(:,:,:,iaa), 'T', 1._wp, kfillmode = jpfillcopy ) 373 ELSE ; CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fnow(:,:,: ), & 374 & sdjf%nrec(1,iaa), sdjf%lsmname ) 375 CALL lbc_lnk( 'fldread', sdjf%fnow(:,:,: ), 'T', 1._wp, kfillmode = jpfillcopy ) 368 376 ENDIF 369 377 ELSE … … 382 390 IF( lk_c1d .AND. lmoor ) THEN 383 391 IF( sdjf%ln_tint ) THEN 384 CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fdta(2,2,1, 2), sdjf%nrec_a(1) )385 CALL lbc_lnk( 'fldread', sdjf%fdta(:,:,1, 2),'Z',1.)392 CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fdta(2,2,1,iaa), sdjf%nrec(1,iaa) ) 393 CALL lbc_lnk( 'fldread', sdjf%fdta(:,:,1,iaa),'T',1., kfillmode = jpfillcopy ) 386 394 ELSE 387 CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fnow(2,2,1 ), sdjf%nrec_a(1) )388 CALL lbc_lnk( 'fldread', sdjf%fnow(:,:,1 ),'Z',1.)395 CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fnow(2,2,1 ), sdjf%nrec(1,iaa) ) 396 CALL lbc_lnk( 'fldread', sdjf%fnow(:,:,1 ),'T',1., kfillmode = jpfillcopy ) 389 397 ENDIF 390 398 ELSE 391 IF( sdjf%ln_tint ) THEN ; CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,1,2), sdjf%nrec_a(1) ) 392 ELSE ; CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,1 ), sdjf%nrec_a(1) ) 399 IF( sdjf%ln_tint ) THEN 400 CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,1,iaa), sdjf%nrec(1,iaa), kfill = jpfillcopy ) 401 ELSE 402 CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,1 ), sdjf%nrec(1,iaa), kfill = jpfillcopy ) 393 403 ENDIF 394 404 ENDIF … … 396 406 IF(lk_c1d .AND. lmoor ) THEN 397 407 IF( sdjf%ln_tint ) THEN 398 CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fdta(2,2,:, 2), sdjf%nrec_a(1) )399 CALL lbc_lnk( 'fldread', sdjf%fdta(:,:,:, 2),'Z',1.)408 CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fdta(2,2,:,iaa), sdjf%nrec(1,iaa) ) 409 CALL lbc_lnk( 'fldread', sdjf%fdta(:,:,:,iaa),'T',1., kfillmode = jpfillcopy ) 400 410 ELSE 401 CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fnow(2,2,: ), sdjf%nrec_a(1) )402 CALL lbc_lnk( 'fldread', sdjf%fnow(:,:,: ),'Z',1.)411 CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fnow(2,2,: ), sdjf%nrec(1,iaa) ) 412 CALL lbc_lnk( 'fldread', sdjf%fnow(:,:,: ),'T',1., kfillmode = jpfillcopy ) 403 413 ENDIF 404 414 ELSE 405 IF( sdjf%ln_tint ) THEN ; CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) ) 406 ELSE ; CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1) ) 415 IF( sdjf%ln_tint ) THEN 416 CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,:,iaa), sdjf%nrec(1,iaa), kfill = jpfillcopy ) 417 ELSE 418 CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec(1,iaa), kfill = jpfillcopy ) 407 419 ENDIF 408 420 ENDIF … … 410 422 ENDIF 411 423 ! 412 sdjf%rotn( 2) = .false. ! vector not yet rotated424 sdjf%rotn(iaa) = .false. ! vector not yet rotated 413 425 ! 414 426 END SUBROUTINE fld_get … … 941 953 IF( sdjf%num <= 0 .OR. .NOT. sdjf%ln_clim ) THEN 942 954 IF( sdjf%num > 0 ) CALL iom_close( sdjf%num ) ! close file if already open 943 CALL iom_open( sdjf%clname, sdjf%num, ldstop = llstop, ldiof = LEN (TRIM(sdjf%wgtname)) > 0 )955 CALL iom_open( sdjf%clname, sdjf%num, ldstop = llstop, ldiof = LEN_TRIM(sdjf%wgtname) > 0 ) 944 956 ENDIF 945 957 ! … … 963 975 ENDIF 964 976 ! 965 CALL iom_open( sdjf%clname, sdjf%num, ldiof = LEN (TRIM(sdjf%wgtname)) > 0 )977 CALL iom_open( sdjf%clname, sdjf%num, ldiof = LEN_TRIM(sdjf%wgtname) > 0 ) 966 978 ! 967 979 ENDIF … … 998 1010 sdf(jf)%cltype = sdf_n(jf)%cltype 999 1011 sdf(jf)%num = -1 1012 sdf(jf)%nbb = 1 ! start with before data in 1 1013 sdf(jf)%naa = 2 ! start with after data in 2 1000 1014 sdf(jf)%wgtname = " " 1001 1015 IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) sdf(jf)%wgtname = TRIM( cdir )//sdf_n(jf)%wname … … 1050 1064 !!---------------------------------------------------------------------- 1051 1065 TYPE( FLD ), INTENT(in ) :: sd ! field with name of weights file 1052 INTEGER , INTENT( inout) :: kwgt ! index of weights1066 INTEGER , INTENT( out) :: kwgt ! index of weights 1053 1067 ! 1054 1068 INTEGER :: kw, nestid ! local integer 1055 LOGICAL :: found ! local logical1056 1069 !!---------------------------------------------------------------------- 1057 1070 ! 1058 1071 !! search down linked list 1059 1072 !! weights filename is either present or we hit the end of the list 1060 found = .FALSE.1061 1073 ! 1062 1074 !! because agrif nest part of filenames are now added in iom_open … … 1068 1080 #endif 1069 1081 DO kw = 1, nxt_wgt-1 1070 IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname).AND. &1071 ref_wgts(kw)%nestid == nestid) THEN1082 IF( ref_wgts(kw)%wgtname == sd%wgtname .AND. & 1083 ref_wgts(kw)%nestid == nestid) THEN 1072 1084 kwgt = kw 1073 found = .TRUE. 1074 EXIT 1085 RETURN 1075 1086 ENDIF 1076 1087 END DO 1077 IF( .NOT.found ) THEN 1078 kwgt = nxt_wgt 1079 CALL fld_weight( sd ) 1080 ENDIF 1088 kwgt = nxt_wgt 1089 CALL fld_weight( sd ) 1081 1090 ! 1082 1091 END SUBROUTINE wgt_list … … 1121 1130 TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file 1122 1131 !! 1123 INTEGER :: j n! dummy loop indices1132 INTEGER :: ji,jj,jn ! dummy loop indices 1124 1133 INTEGER :: inum ! local logical unit 1125 1134 INTEGER :: id ! local variable id … … 1127 1136 INTEGER :: zwrap ! local integer 1128 1137 LOGICAL :: cyclical ! 1129 CHARACTER (len=5) :: aname !1130 INTEGER , DIMENSION( :), ALLOCATABLE:: ddims1131 INTEGER , DIMENSION(jpi,jpj) :: data_src1138 CHARACTER (len=5) :: clname ! 1139 INTEGER , DIMENSION(4) :: ddims 1140 INTEGER :: isrc 1132 1141 REAL(wp), DIMENSION(jpi,jpj) :: data_tmp 1133 1142 !!---------------------------------------------------------------------- … … 1142 1151 !! current weights file 1143 1152 1144 !! open input data file (non-model grid) 1145 CALL iom_open( sd%clname, inum, ldiof = LEN(TRIM(sd%wgtname)) > 0 ) 1146 1147 !! get dimensions: we consider 2D data as 3D data with vertical dim size = 1 1148 IF( SIZE(sd%fnow, 3) > 0 ) THEN 1149 ALLOCATE( ddims(4) ) 1150 ELSE 1151 ALLOCATE( ddims(3) ) 1152 ENDIF 1153 id = iom_varid( inum, sd%clvar, ddims ) 1154 1155 !! close it 1156 CALL iom_close( inum ) 1153 !! get data grid dimensions 1154 id = iom_varid( sd%num, sd%clvar, ddims ) 1157 1155 1158 1156 !! now open the weights file 1159 1160 1157 CALL iom_open ( sd%wgtname, inum ) ! interpolation weights 1161 1158 IF( inum > 0 ) THEN … … 1193 1190 !! two possible cases: bilinear (4 weights) or bicubic (16 weights) 1194 1191 id = iom_varid(inum, 'src05', ldstop=.FALSE.) 1195 IF( id <= 0) THEN 1196 ref_wgts(nxt_wgt)%numwgt = 4 1197 ELSE 1198 ref_wgts(nxt_wgt)%numwgt = 16 1199 ENDIF 1200 1201 ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) ) 1202 ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) ) 1203 ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) ) 1192 IF( id <= 0 ) THEN ; ref_wgts(nxt_wgt)%numwgt = 4 1193 ELSE ; ref_wgts(nxt_wgt)%numwgt = 16 1194 ENDIF 1195 1196 ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(Nis0:Nie0,Njs0:Nje0,4) ) 1197 ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(Nis0:Nie0,Njs0:Nje0,4) ) 1198 ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(Nis0:Nie0,Njs0:Nje0,ref_wgts(nxt_wgt)%numwgt) ) 1204 1199 1205 1200 DO jn = 1,4 1206 aname = ' ' 1207 WRITE(aname,'(a3,i2.2)') 'src',jn 1208 data_tmp(:,:) = 0 1209 CALL iom_get ( inum, jpdom_global, aname, data_tmp(:,:) ) 1210 data_src(:,:) = INT(data_tmp(:,:)) 1211 ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1) 1212 ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1) 1201 WRITE(clname,'(a3,i2.2)') 'src',jn 1202 CALL iom_get ( inum, jpdom_global, clname, data_tmp(:,:), cd_type = 'Z' ) ! no call to lbc_lnk 1203 DO_2D_00_00 1204 !!$ isrc = NINT(data_tmp(ji,jj)) - 1 1205 isrc = INT(data_tmp(ji,jj)) - 1 1206 ref_wgts(nxt_wgt)%data_jpi(ji,jj,jn) = 1 + MOD(isrc, ref_wgts(nxt_wgt)%ddims(1)) 1207 ref_wgts(nxt_wgt)%data_jpj(ji,jj,jn) = 1 + isrc / ref_wgts(nxt_wgt)%ddims(1) 1208 END_2D 1213 1209 END DO 1214 1210 1215 1211 DO jn = 1, ref_wgts(nxt_wgt)%numwgt 1216 aname = ' ' 1217 WRITE(aname,'(a3,i2.2)') 'wgt',jn 1218 ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0 1219 CALL iom_get ( inum, jpdom_global, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) ) 1212 WRITE(clname,'(a3,i2.2)') 'wgt',jn 1213 CALL iom_get ( inum, jpdom_global, clname, data_tmp(:,:), cd_type = 'Z' ) ! no call to lbc_lnk 1214 DO_2D_00_00 1215 ref_wgts(nxt_wgt)%data_wgt(ji,jj,jn) = data_tmp(ji,jj) 1216 END_2D 1220 1217 END DO 1221 1218 CALL iom_close (inum) 1222 1219 1223 1220 ! find min and max indices in grid 1224 ref_wgts(nxt_wgt)%botleft( 1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))1225 ref_wgts(nxt_wgt)%botleft( 2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))1221 ref_wgts(nxt_wgt)%botleft( 1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:)) 1222 ref_wgts(nxt_wgt)%botleft( 2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:)) 1226 1223 ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:)) 1227 1224 ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:)) … … 1247 1244 CALL ctl_stop( ' fld_weight : unable to read the file ' ) 1248 1245 ENDIF 1249 1250 DEALLOCATE (ddims )1251 1246 ! 1252 1247 END SUBROUTINE fld_weight … … 1358 1353 1359 1354 1360 SUBROUTINE fld_interp( num, clvar, kw, kk, dta, & 1361 & nrec, lsmfile) 1355 SUBROUTINE fld_interp( num, clvar, kw, kk, dta, nrec, lsmfile) 1362 1356 !!--------------------------------------------------------------------- 1363 1357 !! *** ROUTINE fld_interp *** … … 1377 1371 INTEGER, DIMENSION(3) :: rec1_lsm, recn_lsm ! temporary arrays for start and length in case of seaoverland 1378 1372 INTEGER :: ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2 ! temporary indices 1379 INTEGER :: jk, jn, jm, jir, jjr ! loop counters 1373 INTEGER :: ji, jj, jk, jn, jir, jjr ! loop counters 1374 INTEGER :: ipk 1380 1375 INTEGER :: ni, nj ! lengths 1381 1376 INTEGER :: jpimin,jpiwid ! temporary indices … … 1388 1383 REAL(wp),DIMENSION(:,:,:), ALLOCATABLE :: ztmp_fly_dta ! local array of values on input grid 1389 1384 !!---------------------------------------------------------------------- 1385 ipk = SIZE(dta, 3) 1390 1386 ! 1391 1387 !! for weighted interpolation we have weights at four corners of a box surrounding … … 1417 1413 1418 1414 1419 IF( LEN ( TRIM(lsmfile)) > 0 ) THEN1415 IF( LEN_TRIM(lsmfile) > 0 ) THEN 1420 1416 !! indeces for ztmp_fly_dta 1421 1417 ! -------------------------- … … 1476 1472 !! first four weights common to both bilinear and bicubic 1477 1473 !! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft 1478 !! note that we have to offset by 1 into fly_dta array because of halo 1479 dta(:,:,:) = 0.0 1480 DO jk = 1,4 1481 DO jn = 1, jpj 1482 DO jm = 1,jpi 1483 ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 1484 nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 1485 dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,:) 1486 END DO 1487 END DO 1474 !! note that we have to offset by 1 into fly_dta array because of halo added to fly_dta (rec1 definition) 1475 dta(:,:,:) = 0._wp 1476 DO jn = 1,4 1477 DO_3D_00_00( 1,ipk ) 1478 ni = ref_wgts(kw)%data_jpi(ji,jj,jn) + 1 1479 nj = ref_wgts(kw)%data_jpj(ji,jj,jn) + 1 1480 dta(ji,jj,jk) = dta(ji,jj,jk) + ref_wgts(kw)%data_wgt(ji,jj,jn) * ref_wgts(kw)%fly_dta(ni,nj,jk) 1481 END_3D 1488 1482 END DO 1489 1483 1490 1484 IF(ref_wgts(kw)%numwgt .EQ. 16) THEN 1491 1485 1492 !! fix up halo points that we couldnt read from file 1493 IF( jpi1 == 2 ) THEN 1494 ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:) 1495 ENDIF 1496 IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 1497 ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:) 1498 ENDIF 1499 IF( jpj1 == 2 ) THEN 1500 ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:) 1501 ENDIF 1502 IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN 1503 ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:) 1504 ENDIF 1505 1506 !! if data grid is cyclic we can do better on east-west edges 1507 !! but have to allow for whether first and last columns are coincident 1508 IF( ref_wgts(kw)%cyclic ) THEN 1509 rec1(2) = MAX( jpjmin-1, 1 ) 1510 recn(1) = 1 1511 recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 ) 1512 jpj1 = 2 + rec1(2) - jpjmin 1513 jpj2 = jpj1 + recn(2) - 1 1514 IF( jpi1 == 2 ) THEN 1515 rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap 1516 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, kstart = rec1, kcount = recn) 1517 ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:) 1518 ENDIF 1519 IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 1520 rec1(1) = 1 + ref_wgts(kw)%overlap 1521 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, kstart = rec1, kcount = recn) 1522 ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:) 1523 ENDIF 1524 ENDIF 1525 1526 ! gradient in the i direction 1527 DO jk = 1,4 1528 DO jn = 1, jpj 1529 DO jm = 1,jpi 1530 ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 1531 nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 1532 dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 * & 1533 (ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:)) 1534 END DO 1535 END DO 1536 END DO 1537 1538 ! gradient in the j direction 1539 DO jk = 1,4 1540 DO jn = 1, jpj 1541 DO jm = 1,jpi 1542 ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 1543 nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 1544 dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 * & 1545 (ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:)) 1546 END DO 1547 END DO 1548 END DO 1549 1550 ! gradient in the ij direction 1551 DO jk = 1,4 1552 DO jn = 1, jpj 1553 DO jm = 1,jpi 1554 ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 1555 nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 1556 dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( & 1557 (ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni ,nj+2,:)) - & 1558 (ref_wgts(kw)%fly_dta(ni+2,nj ,:) - ref_wgts(kw)%fly_dta(ni ,nj ,:))) 1559 END DO 1560 END DO 1486 !! fix up halo points that we couldnt read from file 1487 IF( jpi1 == 2 ) THEN 1488 ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:) 1489 ENDIF 1490 IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 1491 ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:) 1492 ENDIF 1493 IF( jpj1 == 2 ) THEN 1494 ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:) 1495 ENDIF 1496 IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .LT. jpjwid+2 ) THEN 1497 ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:) 1498 ENDIF 1499 1500 !! if data grid is cyclic we can do better on east-west edges 1501 !! but have to allow for whether first and last columns are coincident 1502 IF( ref_wgts(kw)%cyclic ) THEN 1503 rec1(2) = MAX( jpjmin-1, 1 ) 1504 recn(1) = 1 1505 recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 ) 1506 jpj1 = 2 + rec1(2) - jpjmin 1507 jpj2 = jpj1 + recn(2) - 1 1508 IF( jpi1 == 2 ) THEN 1509 rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap 1510 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, kstart = rec1, kcount = recn) 1511 ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:) 1512 ENDIF 1513 IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 1514 rec1(1) = 1 + ref_wgts(kw)%overlap 1515 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, kstart = rec1, kcount = recn) 1516 ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:) 1517 ENDIF 1518 ENDIF 1519 ! 1520 !!$ DO jn = 1,4 1521 !!$ DO_3D_00_00( 1,ipk ) 1522 !!$ ni = ref_wgts(kw)%data_jpi(ji,jj,jn) + 1 1523 !!$ nj = ref_wgts(kw)%data_jpj(ji,jj,jn) + 1 1524 !!$ dta(ji,jj,jk) = dta(ji,jj,jk) & 1525 !!$ ! gradient in the i direction 1526 !!$ & + ref_wgts(kw)%data_wgt(ji,jj,jn+4) * 0.5_wp * & 1527 !!$ & (ref_wgts(kw)%fly_dta(ni+1,nj ,jk) - ref_wgts(kw)%fly_dta(ni-1,nj ,jk)) & 1528 !!$ ! gradient in the j direction 1529 !!$ & + ref_wgts(kw)%data_wgt(ji,jj,jn+8) * 0.5_wp * & 1530 !!$ & (ref_wgts(kw)%fly_dta(ni ,nj+1,jk) - ref_wgts(kw)%fly_dta(ni ,nj-1,jk)) & 1531 !!$ ! gradient in the ij direction 1532 !!$ & + ref_wgts(kw)%data_wgt(ji,jj,jn+12) * 0.25_wp * & 1533 !!$ & ((ref_wgts(kw)%fly_dta(ni+1,nj+1,jk) - ref_wgts(kw)%fly_dta(ni-1,nj+1,jk)) - & 1534 !!$ & (ref_wgts(kw)%fly_dta(ni+1,nj-1,jk) - ref_wgts(kw)%fly_dta(ni-1,nj-1,jk))) 1535 !!$ END_3D 1536 !!$ END DO 1537 ! 1538 DO jn = 1,4 1539 DO_3D_00_00( 1,ipk ) 1540 ni = ref_wgts(kw)%data_jpi(ji,jj,jn) 1541 nj = ref_wgts(kw)%data_jpj(ji,jj,jn) 1542 ! gradient in the i direction 1543 dta(ji,jj,jk) = dta(ji,jj,jk) + ref_wgts(kw)%data_wgt(ji,jj,jn+4) * 0.5_wp * & 1544 & (ref_wgts(kw)%fly_dta(ni+2,nj+1,jk) - ref_wgts(kw)%fly_dta(ni ,nj+1,jk)) 1545 END_3D 1546 END DO 1547 DO jn = 1,4 1548 DO_3D_00_00( 1,ipk ) 1549 ni = ref_wgts(kw)%data_jpi(ji,jj,jn) 1550 nj = ref_wgts(kw)%data_jpj(ji,jj,jn) 1551 ! gradient in the j direction 1552 dta(ji,jj,jk) = dta(ji,jj,jk) + ref_wgts(kw)%data_wgt(ji,jj,jn+8) * 0.5_wp * & 1553 & (ref_wgts(kw)%fly_dta(ni+1,nj+2,jk) - ref_wgts(kw)%fly_dta(ni+1,nj ,jk)) 1554 END_3D 1555 END DO 1556 DO jn = 1,4 1557 DO_3D_00_00( 1,ipk ) 1558 ni = ref_wgts(kw)%data_jpi(ji,jj,jn) 1559 nj = ref_wgts(kw)%data_jpj(ji,jj,jn) 1560 ! gradient in the ij direction 1561 dta(ji,jj,jk) = dta(ji,jj,jk) + ref_wgts(kw)%data_wgt(ji,jj,jn+12) * 0.25_wp * ( & 1562 & (ref_wgts(kw)%fly_dta(ni+2,nj+2,jk) - ref_wgts(kw)%fly_dta(ni ,nj+2,jk)) - & 1563 & (ref_wgts(kw)%fly_dta(ni+2,nj ,jk) - ref_wgts(kw)%fly_dta(ni ,nj ,jk))) 1564 END_3D 1561 1565 END DO 1562 1566 !
Note: See TracChangeset
for help on using the changeset viewer.
