- Timestamp:
- 2013-11-18T12:57:11+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_LOCEAN_CMCC_INGV_2013/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
r3851 r4230 7 7 !! ! 05-2008 (S. Alderson) Modified for Interpolation in memory 8 8 !! ! from input grid to model grid 9 !! ! 10-2013 (D. Delrosso, P. Oddo) implement suppression of 10 !! ! land point prior to interpolation 9 11 !!---------------------------------------------------------------------- 10 12 … … 22 24 USE wrk_nemo ! work arrays 23 25 USE ioipsl, ONLY : ymds2ju, ju2ymds ! for calendar 24 26 USE sbc_oce 27 25 28 IMPLICIT NONE 26 29 PRIVATE … … 40 43 ! ! a string starting with "U" or "V" for each component 41 44 ! ! chars 2 onwards identify which components go together 45 CHARACTER(len = 34) :: lname ! generic name of a NetCDF land/sea mask file to be used, blank if not 46 ! ! 0=sea 1=land 42 47 END TYPE FLD_N 43 48 … … 60 65 LOGICAL, DIMENSION(2) :: rotn ! flag to indicate whether before/after field has been rotated 61 66 INTEGER :: nreclast ! last record to be read in the current file 67 CHARACTER(len = 256) :: lsmname ! current name of the NetCDF mask file acting as a key 62 68 END TYPE FLD 63 69 … … 95 101 TYPE( WGT ), DIMENSION(tot_wgts) :: ref_wgts ! array of wgts 96 102 INTEGER :: nxt_wgt = 1 ! point to next available space in ref_wgts array 103 REAL(wp), PARAMETER :: undeff_lsm = -999.00_wp 97 104 98 105 !$AGRIF_END_DO_NOT_TREAT … … 591 598 ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 592 599 CALL wgt_list( sdjf, iw ) 593 IF( sdjf%ln_tint ) THEN ; CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) ) 594 ELSE ; CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fnow(:,:,: ), sdjf%nrec_a(1) ) 600 IF( sdjf%ln_tint ) THEN ; CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fdta(:,:,:,2), & 601 & sdjf%nrec_a(1), sdjf%lsmname ) 602 ELSE ; CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fnow(:,:,: ), & 603 & sdjf%nrec_a(1), sdjf%lsmname ) 595 604 ENDIF 596 605 ELSE … … 856 865 sdf(jf)%wgtname = " " 857 866 IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname ) 867 sdf(jf)%lsmname = " " 868 IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 ) sdf(jf)%lsmname = TRIM( cdir )//TRIM( sdf_n(jf)%lname ) 858 869 sdf(jf)%vcomp = sdf_n(jf)%vcomp 859 870 sdf(jf)%rotn(:) = .TRUE. ! pretend to be rotated -> won't try to rotate data before the first call to fld_get … … 878 889 & ' weights : ' , TRIM( sdf(jf)%wgtname ), & 879 890 & ' pairing : ' , TRIM( sdf(jf)%vcomp ), & 880 & ' data type: ' , sdf(jf)%cltype 891 & ' data type: ' , sdf(jf)%cltype , & 892 & ' land/sea mask:' , TRIM( sdf(jf)%lsmname ) 881 893 call flush(numout) 882 894 END DO … … 1098 1110 1099 1111 1100 SUBROUTINE fld_interp( num, clvar, kw, kk, dta, nrec ) 1112 SUBROUTINE apply_seaoverland(clmaskfile,zfieldo,jpi1_lsm,jpi2_lsm,jpj1_lsm, & 1113 & jpj2_lsm,itmpi,itmpj,itmpz,rec1_lsm,recn_lsm) 1114 !!--------------------------------------------------------------------- 1115 !! *** ROUTINE apply_seaoverland *** 1116 !! 1117 !! ** Purpose : avoid spurious fluxes in coastal or near-coastal areas 1118 !! due to the wrong usage of "land" values from the coarse 1119 !! atmospheric model when spatial interpolation is required 1120 !! D. Delrosso INGV 1121 !!---------------------------------------------------------------------- 1122 INTEGER :: inum,jni,jnj,jnz,jc ! temporary indices 1123 INTEGER, INTENT(in) :: itmpi,itmpj,itmpz ! lengths 1124 INTEGER, INTENT(in) :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices 1125 INTEGER, DIMENSION(3), INTENT(in) :: rec1_lsm,recn_lsm ! temporary arrays for start and length 1126 REAL(wp),DIMENSION (:,:,:),INTENT(inout) :: zfieldo ! input/output array for seaoverland application 1127 REAL(wp),DIMENSION (:,:,:),ALLOCATABLE :: zslmec1 ! temporary array for land point detection 1128 REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfieldn ! array of forcing field with undeff for land points 1129 REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfield ! array of forcing field 1130 CHARACTER (len=100), INTENT(in) :: clmaskfile ! land/sea mask file name 1131 !!--------------------------------------------------------------------- 1132 ALLOCATE ( zslmec1(itmpi,itmpj,itmpz) ) 1133 ALLOCATE ( zfieldn(itmpi,itmpj) ) 1134 ALLOCATE ( zfield(itmpi,itmpj) ) 1135 1136 ! Retrieve the land sea mask data 1137 CALL iom_open( clmaskfile, inum ) 1138 SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) ) 1139 CASE(1) 1140 CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm) 1141 CASE DEFAULT 1142 CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm) 1143 END SELECT 1144 CALL iom_close( inum ) 1145 1146 DO jnz=1,rec1_lsm(3) !! Loop over k dimension 1147 1148 DO jni=1,itmpi !! copy the original field into a tmp array 1149 DO jnj=1,itmpj !! substituting undeff over land points 1150 zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz) 1151 IF ( zslmec1(jni,jnj,jnz) == 1. ) THEN 1152 zfieldn(jni,jnj) = undeff_lsm 1153 ENDIF 1154 END DO 1155 END DO 1156 1157 CALL seaoverland(zfieldn,itmpi,itmpj,zfield) 1158 DO jc=1,nn_lsm 1159 CALL seaoverland(zfield,itmpi,itmpj,zfield) 1160 END DO 1161 1162 ! Check for Undeff and substitute original values 1163 IF(ANY(zfield==undeff_lsm)) THEN 1164 DO jni=1,itmpi 1165 DO jnj=1,itmpj 1166 IF (zfield(jni,jnj)==undeff_lsm) THEN 1167 zfield(jni,jnj) = zfieldo(jni,jnj,jnz) 1168 ENDIF 1169 ENDDO 1170 ENDDO 1171 ENDIF 1172 1173 zfieldo(:,:,jnz)=zfield(:,:) 1174 1175 END DO !! End Loop over k dimension 1176 1177 DEALLOCATE ( zslmec1 ) 1178 DEALLOCATE ( zfieldn ) 1179 DEALLOCATE ( zfield ) 1180 1181 END SUBROUTINE apply_seaoverland 1182 1183 1184 SUBROUTINE seaoverland(zfieldn,ileni,ilenj,zfield) 1185 !!--------------------------------------------------------------------- 1186 !! *** ROUTINE seaoverland *** 1187 !! 1188 !! ** Purpose : create shifted matrices for seaoverland application 1189 !! D. Delrosso INGV 1190 !!---------------------------------------------------------------------- 1191 INTEGER,INTENT(in) :: ileni,ilenj ! lengths 1192 REAL,DIMENSION (ileni,ilenj),INTENT(in) :: zfieldn ! array of forcing field with undeff for land points 1193 REAL,DIMENSION (ileni,ilenj),INTENT(out) :: zfield ! array of forcing field 1194 REAL,DIMENSION (ileni,ilenj) :: zmat1,zmat2,zmat3,zmat4 ! temporary arrays for seaoverland application 1195 REAL,DIMENSION (ileni,ilenj) :: zmat5,zmat6,zmat7,zmat8 ! temporary arrays for seaoverland application 1196 REAL,DIMENSION (ileni,ilenj) :: zlsm2d ! temporary arrays for seaoverland application 1197 REAL,DIMENSION (ileni,ilenj,8) :: zlsm3d ! temporary arrays for seaoverland application 1198 LOGICAL,DIMENSION (ileni,ilenj,8) :: ll_msknan3d ! logical mask for undeff detection 1199 LOGICAL,DIMENSION (ileni,ilenj) :: ll_msknan2d ! logical mask for undeff detection 1200 !!---------------------------------------------------------------------- 1201 zmat8 = eoshift(zfieldn , SHIFT=-1, BOUNDARY = (/zfieldn(:,1)/) ,DIM=2) 1202 zmat1 = eoshift(zmat8 , SHIFT=-1, BOUNDARY = (/zmat8(1,:)/) ,DIM=1) 1203 zmat2 = eoshift(zfieldn , SHIFT=-1, BOUNDARY = (/zfieldn(1,:)/) ,DIM=1) 1204 zmat4 = eoshift(zfieldn , SHIFT= 1, BOUNDARY = (/zfieldn(:,ilenj)/),DIM=2) 1205 zmat3 = eoshift(zmat4 , SHIFT=-1, BOUNDARY = (/zmat4(1,:)/) ,DIM=1) 1206 zmat5 = eoshift(zmat4 , SHIFT= 1, BOUNDARY = (/zmat4(ileni,:)/) ,DIM=1) 1207 zmat6 = eoshift(zfieldn , SHIFT= 1, BOUNDARY = (/zfieldn(ileni,:)/),DIM=1) 1208 zmat7 = eoshift(zmat8 , SHIFT= 1, BOUNDARY = (/zmat8(ileni,:)/) ,DIM=1) 1209 1210 zlsm3d = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /)) 1211 ll_msknan3d = .not.(zlsm3d==undeff_lsm) 1212 ll_msknan2d = .not.(zfieldn==undeff_lsm) ! FALSE where is Undeff (land) 1213 zlsm2d = (SUM ( zlsm3d, 3 , ll_msknan3d ) )/(MAX(1,(COUNT( ll_msknan3d , 3 )) )) 1214 WHERE ((COUNT( ll_msknan3d , 3 )) == 0.0_wp) zlsm2d = undeff_lsm 1215 zfield = MERGE (zfieldn,zlsm2d,ll_msknan2d) 1216 END SUBROUTINE seaoverland 1217 1218 1219 SUBROUTINE fld_interp( num, clvar, kw, kk, dta, & 1220 & nrec, lsmfile) 1101 1221 !!--------------------------------------------------------------------- 1102 1222 !! *** ROUTINE fld_interp *** … … 1111 1231 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: dta ! output field on model grid 1112 1232 INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice) 1233 CHARACTER(LEN=*) , INTENT(in ) :: lsmfile ! land sea mask file name 1113 1234 !! 1114 INTEGER, DIMENSION(3) :: rec1,recn ! temporary arrays for start and length 1115 INTEGER :: jk, jn, jm ! loop counters 1116 INTEGER :: ni, nj ! lengths 1117 INTEGER :: jpimin,jpiwid ! temporary indices 1118 INTEGER :: jpjmin,jpjwid ! temporary indices 1119 INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices 1235 REAL(wp),DIMENSION(:,:,:),ALLOCATABLE :: ztmp_fly_dta,zfieldo ! temporary array of values on input grid 1236 INTEGER, DIMENSION(3) :: rec1,recn ! temporary arrays for start and length 1237 INTEGER, DIMENSION(3) :: rec1_lsm,recn_lsm ! temporary arrays for start and length in case of seaoverland 1238 INTEGER :: ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2 ! temporary indices 1239 INTEGER :: jk, jn, jm, jir, jjr ! loop counters 1240 INTEGER :: ni, nj ! lengths 1241 INTEGER :: jpimin,jpiwid ! temporary indices 1242 INTEGER :: jpimin_lsm,jpiwid_lsm ! temporary indices 1243 INTEGER :: jpjmin,jpjwid ! temporary indices 1244 INTEGER :: jpjmin_lsm,jpjwid_lsm ! temporary indices 1245 INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices 1246 INTEGER :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices 1247 INTEGER :: itmpi,itmpj,itmpz ! lengths 1248 1120 1249 !!---------------------------------------------------------------------- 1121 1250 ! … … 1147 1276 jpj2 = jpj1 + recn(2) - 1 1148 1277 1149 ref_wgts(kw)%fly_dta(:,:,:) = 0.0 1150 SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) ) 1151 CASE(1) 1152 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn) 1153 CASE DEFAULT 1154 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn) 1155 END SELECT 1278 1279 IF( LEN( TRIM(lsmfile) ) > 0 ) THEN 1280 !! indeces for ztmp_fly_dta 1281 ! -------------------------- 1282 rec1_lsm(1)=MAX(rec1(1)-nn_lsm,1) ! starting index for enlarged external data, x direction 1283 rec1_lsm(2)=MAX(rec1(2)-nn_lsm,1) ! starting index for enlarged external data, y direction 1284 rec1_lsm(3) = 1 ! vertical dimension 1285 recn_lsm(1)=MIN(rec1(1)-rec1_lsm(1)+recn(1)+nn_lsm,ref_wgts(kw)%ddims(1)-rec1_lsm(1)) ! n points in x direction 1286 recn_lsm(2)=MIN(rec1(2)-rec1_lsm(2)+recn(2)+nn_lsm,ref_wgts(kw)%ddims(2)-rec1_lsm(2)) ! n points in y direction 1287 recn_lsm(3) = kk ! number of vertical levels in the input file 1288 1289 ! Avoid out of bound 1290 jpimin_lsm = MAX( rec1_lsm(1)+1, 1 ) 1291 jpjmin_lsm = MAX( rec1_lsm(2)+1, 1 ) 1292 jpiwid_lsm = MIN( recn_lsm(1)-2,ref_wgts(kw)%ddims(1)-rec1(1)+1) 1293 jpjwid_lsm = MIN( recn_lsm(2)-2,ref_wgts(kw)%ddims(2)-rec1(2)+1) 1294 1295 jpi1_lsm = 2+rec1_lsm(1)-jpimin_lsm 1296 jpj1_lsm = 2+rec1_lsm(2)-jpjmin_lsm 1297 jpi2_lsm = jpi1_lsm + recn_lsm(1) - 1 1298 jpj2_lsm = jpj1_lsm + recn_lsm(2) - 1 1299 1300 1301 itmpi=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),1) 1302 itmpj=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),2) 1303 itmpz=kk 1304 ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz)) 1305 ztmp_fly_dta(:,:,:) = 0.0 1306 SELECT CASE( SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) ) 1307 CASE(1) 1308 CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), & 1309 & nrec, rec1_lsm, recn_lsm) 1310 CASE DEFAULT 1311 CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), & 1312 & nrec, rec1_lsm, recn_lsm) 1313 END SELECT 1314 CALL apply_seaoverland(lsmfile,ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), & 1315 & jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm, & 1316 & itmpi,itmpj,itmpz,rec1_lsm,recn_lsm) 1317 1318 1319 ! Relative indeces for remapping 1320 ii_lsm1 = (rec1(1)-rec1_lsm(1))+1 1321 ii_lsm2 = (ii_lsm1+recn(1))-1 1322 ij_lsm1 = (rec1(2)-rec1_lsm(2))+1 1323 ij_lsm2 = (ij_lsm1+recn(2))-1 1324 1325 ref_wgts(kw)%fly_dta(:,:,:) = 0.0 1326 ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:) = ztmp_fly_dta(ii_lsm1:ii_lsm2,ij_lsm1:ij_lsm2,:) 1327 DEALLOCATE(ztmp_fly_dta) 1328 1329 ELSE 1330 1331 ref_wgts(kw)%fly_dta(:,:,:) = 0.0 1332 SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) ) 1333 CASE(1) 1334 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn) 1335 CASE DEFAULT 1336 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn) 1337 END SELECT 1338 ENDIF 1339 1156 1340 1157 1341 !! first four weights common to both bilinear and bicubic
Note: See TracChangeset
for help on using the changeset viewer.
