|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 SUBROUTINE ASSIM_NATURE_ISBA_OI(YPROGRAM, KI, & 00002 PRRCL, PRRSL, PRRCN, PRRSN, & 00003 PATMNEB, PITM, PEVAPTR, PEVAP, & 00004 PSNC, PTSC, & 00005 PTS_O, PT2M_O, PHU2M_O, PSWE, & 00006 HTEST ) 00007 00008 ! ------------------------------------------------------------------------------------------ 00009 ! ***************************************************************************************** 00010 ! 00011 ! Routine to perform OI within SURFEX 00012 ! a soil analysis for water content and temperature 00013 ! using the Meteo-France optimum interpolation technique of Giard and Bazile (2000) 00014 ! 00015 ! Derived from CANARI subroutines externalized by Lora Taseva (Dec. 2007) 00016 ! 00017 ! Author : Jean-Francois Mahfouf (01/2008) 00018 ! 00019 ! Modifications : 00020 ! (05/2008) : The I/O of this version follow the newly available LFI format in SURFEX 00021 ! (01/2009) : Read directly atmospheric FA files using XRD library instead of using "edf" 00022 ! (06/2009) : Modifications to allow the assimilation of ASCAT superficial soil moisture 00023 ! (09/2010) : More parameters to goto_surfex 00024 ! (03/2011) : Initialization of ZEVAPTR (F.Bouyssel) 00025 ! (07/2011) : Read pgd+prep (B. Decharme) 00026 ! (04/2012) : Made as a subroutine (T. Aspelien) 00027 ! ****************************************************************************************** 00028 ! ------------------------------------------------------------------------------------------ 00029 USE MODD_TYPE_DATE_SURF, ONLY : DATE_TIME 00030 USE MODD_CSTS, ONLY : XDAY, XPI, XG, XRHOLW, XLVTT, NDAYSEC 00031 USE MODD_SURF_PAR, ONLY : XUNDEF 00032 ! 00033 USE MODD_ASSIM, ONLY : LOBSWG,ITRAD,LPRINT,NECHGU,RD1,RSCALDW,& 00034 RTHR_QC,SIGWGB,SIGWGO,SIGWGO_MAX 00035 USE MODN_IO_OFFLINE, ONLY : CPGDFILE,CPREPFILE 00036 ! 00037 #ifdef LFI 00038 USE MODD_IO_SURF_LFI, ONLY : CFILEIN_LFI, CFILEOUT_LFI 00039 #endif 00040 ! 00041 USE MODD_SURF_ATM_n, ONLY : NDIM_NATURE,NR_NATURE 00042 USE MODD_SURF_ATM_GRID_n, ONLY : XLAT, XLON 00043 ! 00044 USE YOMHOOK, ONLY : LHOOK, DR_HOOK 00045 USE PARKIND1, ONLY : JPRB 00046 ! 00047 USE MODI_ABOR1_SFX 00048 USE MODI_INIT_IO_SURF_n 00049 USE MODI_READ_SURF 00050 USE MODI_END_IO_SURF_n 00051 USE MODI_TRANS_CHAINE 00052 USE MODI_IO_BUFF_CLEAN_n 00053 USE MODI_OI_BC_SOIL_MOISTURE 00054 USE MODI_OI_CACSTS 00055 USE MODI_FLAG_UPDATE 00056 USE MODI_WRITE_SURF 00057 ! 00058 IMPLICIT NONE 00059 CHARACTER(LEN=6), INTENT(IN) :: YPROGRAM ! program calling surf. schemes 00060 INTEGER, INTENT(IN) :: KI 00061 REAL, DIMENSION(KI), INTENT(IN) :: PRRCL 00062 REAL, DIMENSION(KI), INTENT(IN) :: PRRSL 00063 REAL, DIMENSION(KI), INTENT(IN) :: PRRCN 00064 REAL, DIMENSION(KI), INTENT(IN) :: PRRSN 00065 REAL, DIMENSION(KI), INTENT(IN) :: PATMNEB 00066 REAL, DIMENSION(KI), INTENT(IN) :: PITM 00067 REAL, DIMENSION(KI), INTENT(IN) :: PEVAPTR 00068 REAL, DIMENSION(KI), INTENT(IN) :: PEVAP 00069 REAL, DIMENSION(KI), INTENT(IN) :: PSNC 00070 REAL, DIMENSION(KI), INTENT(IN) :: PTSC 00071 REAL, DIMENSION(KI), INTENT(IN) :: PTS_O 00072 REAL, DIMENSION(KI), INTENT(IN) :: PT2M_O 00073 REAL, DIMENSION(KI), INTENT(IN) :: PHU2M_O 00074 REAL, DIMENSION(KI), INTENT(OUT):: PSWE 00075 CHARACTER(LEN=2), INTENT(IN) :: HTEST ! must be equal to 'OK' 00076 ! 00077 INTEGER :: IDAT 00078 ! 00079 ! Declarations of local variables 00080 ! 00081 CHARACTER(LEN=3) :: YREAD 00082 CHARACTER(LEN=2) :: CMONTH 00083 INTEGER :: IYEAR ! current year (UTC) 00084 INTEGER :: IMONTH ! current month (UTC) 00085 INTEGER :: IDAY ! current day (UTC) 00086 INTEGER :: NSSSSS ! current time since start of the run (s) 00087 INTEGER :: IRESP ! return code 00088 TYPE (DATE_TIME) :: TTIME ! Current date and time 00089 00090 ! 00091 ! Arrays for soil OI analysis 00092 ! 00093 00094 REAL, DIMENSION (KI) :: PWS, PWP, PTS, PTP, PTL, PLAI, PVEG, PRSMIN, PD2, PSAB, PARG, 00095 PTSUN, PZENITH, PAZIMSOL, PTCLS, PHCLS, PRAIN, PSNOW, 00096 PWIND, PSWD, PSWS, PUCLS, PVCLS, PSNS, 00097 ZT2INC, ZH2INC, ZWS, ZWP, ZTL, ZTS, ZTP, ZTCLS, ZHCLS, ZUCLS, 00098 ZVCLS, PSSTC, PWPINC1, PWPINC2, PWPINC3, PT2MBIAS, PH2MBIAS, 00099 PALBF, PEMISF, PZ0F, PIVEG, PZ0H, 00100 PTPC, PWSC, PWPC, ZEVAP, ZEVAPTR, PGELAT, PGELAM, PGEMU, 00101 ZWSINC, ZWPINC, ZSNS, ZTSINC, ZTPINC, ZTLINC, 00102 PSM_O, PSIG_SMO, PLSM_O, PWS_O, ZWGINC , 00103 ZALT 00104 00105 REAL,DIMENSION(KI) :: PLON, PLAT 00106 INTEGER :: I,J 00107 INTEGER :: IVERSION, IBUGFIX 00108 CHARACTER(LEN=10) :: YVAR ! Name of the prognostic variable (in LFI file) 00109 CHARACTER(LEN=100) :: YPREFIX ! Prefix of the prognostic variable (in LFI file) 00110 INTEGER :: INOBS ! number of observations 00111 REAL :: ZTHRES 00112 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00113 ! 00114 ! ---------------------------------------------------------------------------------- 00115 ! 00116 IF (LHOOK) CALL DR_HOOK('ASSIM_NATURE_ISBA_OI',0,ZHOOK_HANDLE) 00117 00118 IF (HTEST/='OK') THEN 00119 CALL ABOR1_SFX('ASSIM_ISBA_n: FATAL ERROR DURING ARGUMENT TRANSFER') 00120 END IF 00121 00122 PRINT *,'--------------------------------------------------------------------------' 00123 PRINT *,'| |' 00124 PRINT *,'| ENTER OI_ASSIM |' 00125 PRINT *,'| |' 00126 PRINT *,'--------------------------------------------------------------------------' 00127 00128 ! 00129 ! Update some constants dependant from NACVEG 00130 ! 00131 ! scaling of soil moisture increments when assimilation window is different 00132 ! from 6 hours 00133 RSCALDW = REAL(NECHGU)/6.0 00134 ! half assimilation window in sec 00135 ITRAD = NECHGU*1800 00136 00137 ! 00138 ! 00139 ! READ PGD FILE 00140 !------------------------------------------------------------ 00141 ! 00142 ! File handling definition 00143 ! 00144 #ifdef LFI 00145 CFILEIN_LFI = CPGDFILE ! input PGD file (surface fields) 00146 #endif 00147 ! 00148 ! Read grid dimension for allocation 00149 ! 00150 CALL INIT_IO_SURF_n(YPROGRAM,'NATURE','SURF ','READ ') 00151 00152 CALL READ_SURF(YPROGRAM,'SAND', PSAB, IRESP) 00153 CALL READ_SURF(YPROGRAM,'CLAY', PARG, IRESP) 00154 CALL READ_SURF(YPROGRAM,'ZS', ZALT, IRESP) 00155 ! 00156 CALL END_IO_SURF_n(YPROGRAM) 00157 ! 00158 !------------------------------------------------------------ 00159 ! READ PREP FILE 00160 !------------------------------------------------------------ 00161 ! 00162 ! File handling definition 00163 ! 00164 #ifdef LFI 00165 CFILEIN_LFI = CPREPFILE ! input PREP file (surface fields) 00166 #endif 00167 ! 00168 ! Read grid dimension for allocation 00169 ! 00170 CALL INIT_IO_SURF_n(YPROGRAM,'NATURE','SURF ','READ ') 00171 ! 00172 ! Read prognostic variables 00173 ! 00174 CALL READ_SURF(YPROGRAM,'VERSION',IVERSION,IRESP) 00175 CALL READ_SURF(YPROGRAM,'BUG',IBUGFIX,IRESP) 00176 00177 CALL READ_SURF(YPROGRAM,'WG1', PWS, IRESP) 00178 CALL READ_SURF(YPROGRAM,'WG2', PWP, IRESP) 00179 CALL READ_SURF(YPROGRAM,'TG1', PTS, IRESP) 00180 CALL READ_SURF(YPROGRAM,'TG2', PTP, IRESP) 00181 CALL READ_SURF(YPROGRAM,'WGI2', PTL, IRESP) 00182 IF (IVERSION<7 .OR. IVERSION==7 .AND. IBUGFIX<3) THEN 00183 CALL READ_SURF(YPROGRAM,'WSNOW_VEG1',PSNS, IRESP) 00184 ELSE 00185 CALL READ_SURF(YPROGRAM,'WSN_VEG1',PSNS, IRESP) 00186 ENDIF 00187 ! 00188 CALL READ_SURF(YPROGRAM,'STORAGETYPE', YREAD, IRESP) 00189 IF (YREAD=='ALL') THEN 00190 CALL READ_SURF(YPROGRAM,'T2M', PTCLS, IRESP) 00191 CALL READ_SURF(YPROGRAM,'HU2M', PHCLS, IRESP) 00192 CALL READ_SURF(YPROGRAM,'ZON10M', PUCLS, IRESP) 00193 CALL READ_SURF(YPROGRAM,'MER10M', PVCLS, IRESP) 00194 ENDIF 00195 ! 00196 ! Read diag of constant surface fields (present in PREP file) 00197 ! 00198 CALL READ_SURF(YPROGRAM,'RSMIN', PRSMIN,IRESP) 00199 CALL READ_SURF(YPROGRAM,'DG2', PD2, IRESP) 00200 CALL READ_SURF(YPROGRAM,'LAI', PLAI, IRESP) 00201 CALL READ_SURF(YPROGRAM,'VEG', PVEG, IRESP) 00202 00203 ! Set PIVEG (SURFIND.VEG.DOMI) since it is not available 00204 PIVEG = 0.0 00205 00206 ! Find current time 00207 ! 00208 CALL READ_SURF(YPROGRAM,'DTCUR',TTIME,IRESP) 00209 ! 00210 ! Time initializations 00211 ! 00212 IYEAR = TTIME%TDATE%YEAR 00213 IMONTH = TTIME%TDATE%MONTH 00214 IDAY = TTIME%TDATE%DAY 00215 NSSSSS = TTIME%TIME 00216 IF (NSSSSS > NDAYSEC) NSSSSS = NSSSSS - NDAYSEC 00217 CALL TRANS_CHAINE(CMONTH,IMONTH,2) 00218 IDAT = IYEAR*10000. + IMONTH*100. + IDAY 00219 00220 ! 00221 ! PRINT 00222 ! 00223 IF (LPRINT) THEN 00224 PRINT *,'value in PREP file => WG1 ',SUM(PWS)/KI 00225 PRINT *,'value in PREP file => WG2 ',SUM(PWP)/KI 00226 PRINT *,'value in PREP file => TG1 ',SUM(PTS)/KI 00227 PRINT *,'value in PREP file => TG2 ',SUM(PTP)/KI 00228 PRINT *,'value in PREP file => WGI2 ',SUM(PTL)/KI 00229 PRINT *,'value in PREP file => WSNOW_VEG1',SUM(PSNS)/KI 00230 PRINT *,'value in PREP file => LAI ',SUM(PLAI)/KI 00231 PRINT *,'value in PREP file => VEG ',SUM(PVEG)/KI 00232 PRINT *,'value in PREP file => RSMIN ',SUM(PRSMIN)/KI 00233 PRINT *,'value in PREP file => DATA_DG2 ',SUM(PD2)/KI 00234 PRINT *,'value in PREP file => SAND ',SUM(PSAB)/KI 00235 PRINT *,'value in PREP file => CLAY ',SUM(PARG)/KI 00236 PRINT *,'value in PREP file => ZS ',SUM(ZALT)/KI 00237 ENDIF 00238 ! 00239 CALL END_IO_SURF_n(YPROGRAM) 00240 CALL IO_BUFF_CLEAN_n 00241 ! 00242 ! 00243 ! Read ASCAT SM observations (in percent) 00244 ! 00245 00246 00247 INOBS = 0 00248 IF (LOBSWG) THEN 00249 OPEN(UNIT=111,FILE='ASCAT_SM.DAT') 00250 DO I=1,KI 00251 READ(111,*) PSM_O(I),PSIG_SMO(I),PLSM_O(I) 00252 IF (PLSM_O(I) < 1.0) PSM_O(I) = 999.0 ! data rejection if not on land 00253 IF (PSIG_SMO(I) > SIGWGO_MAX) PSM_O(I) = 999.0 ! data rejection of error too large 00254 IF (PSM_O(I) /= 999.0) INOBS = INOBS + 1 00255 ENDDO 00256 CLOSE(UNIT=111) 00257 PRINT *,'READ ASCAT SM OK' 00258 ELSE 00259 PSM_O(:) = 999.0 00260 PSIG_SMO(:) = 999.0 00261 PLSM_O(:) = 0.0 00262 ENDIF 00263 PRINT *,' NUMBER OF ASCAT OBSERVATIONS AFTER INITIAL CHECKS :: ',INOBS 00264 INOBS = 0 00265 ! 00266 ! Perform bias correction of SM observations 00267 ! 00268 CALL OI_BC_SOIL_MOISTURE(KI,PSM_O,PSAB,PWS_O) 00269 00270 !Set longitudes/latitudes 00271 DO I=1,NDIM_NATURE 00272 PLON(I)=XLON(NR_NATURE(I)) 00273 PLAT(I)=XLAT(NR_NATURE(I)) 00274 ENDDO 00275 00276 ! 00277 ! Screen-level innovations 00278 ! 00279 ZT2INC(:) = PT2M_O(:) - PTCLS(:) 00280 ZH2INC(:) = PHU2M_O(:) - PHCLS(:) 00281 ! 00282 ! Threshold for background check 00283 ! 00284 ZTHRES=RTHR_QC*SQRT(SIGWGO**2 + SIGWGB**2) 00285 ! 00286 ! Superficial soil moisture innovations in (m3/m3) 00287 ! 00288 DO I=1,KI 00289 IF (PWS_O(I) /= 999.0) THEN 00290 ZWGINC(I) = PWS_O(I) - PWS(I) 00291 IF (ABS(ZWGINC(I)) > ZTHRES) THEN 00292 ZWGINC(I) = 0.0 ! background check 00293 ELSE 00294 INOBS = INOBS + 1 00295 ENDIF 00296 ELSE 00297 ZWGINC(I) = 0.0 00298 ENDIF 00299 ENDDO 00300 PRINT *,' NUMBER OF ASCAT OBSERVATIONS AFTER BACKGROUND CHECK :: ',INOBS 00301 ! 00302 PRINT *,' ' 00303 PRINT *,'Mean T2m increments ',SUM(ZT2INC)/KI 00304 PRINT *,'Mean HU2m increments ',SUM(ZH2INC)/KI 00305 PRINT *,' ' 00306 ! 00307 ! Interface (define arrays and perform unit conversions) 00308 ! 00309 PARG(:) = PARG(:)*100.0 00310 PSAB(:) = PSAB(:)*100.0 00311 ! 00312 ZWS(:) = PWS(:)*RD1*XRHOLW ! conversion of m3/m3 -> mm 00313 ZWP(:) = PWP(:)*PD2(:)*XRHOLW ! conversion of m3/m3 -> mm 00314 ZTL(:) = PTL(:)*PD2(:)*XRHOLW ! conversion of m3/m3 -> mm 00315 ZTCLS(:) = PTCLS(:) 00316 ZHCLS(:) = PHCLS(:) 00317 ZUCLS(:) = PUCLS(:) 00318 ZVCLS(:) = PVCLS(:) 00319 00320 ! SST not used in cacsts 00321 PSSTC(:) = 0. 00322 00323 PWPINC1(:) = XUNDEF 00324 PWPINC2(:) = XUNDEF 00325 PWPINC3(:) = XUNDEF 00326 PT2MBIAS(:) = XUNDEF 00327 PH2MBIAS(:) = XUNDEF 00328 00329 ! 00330 ! Sea-ice surface properties 00331 ! 00332 PALBF(:) = XUNDEF 00333 PEMISF(:) = XUNDEF 00334 PZ0F(:) = XUNDEF 00335 PZ0H(:) = XUNDEF 00336 ! 00337 ! Climatological arrays set to missing values 00338 ! 00339 PWSC(:) = XUNDEF 00340 PWPC(:) = XUNDEF 00341 PTPC(:) = XUNDEF 00342 ! 00343 DO I=1,KI 00344 PGELAT(I) = PLAT(I) 00345 PGELAM(I) = PLON(I) 00346 PGEMU(I) = SIN(PLAT(I)*XPI/180.) 00347 ENDDO 00348 ! 00349 ZEVAP(:) = (PEVAP(:)/XLVTT*XDAY)/(NECHGU*3600.) ! conversion W/m2 -> mm/day 00350 ZEVAPTR(:) = PEVAPTR(:)*XDAY 00351 ZSNS(:) = PSNS(:) 00352 ! 00353 DO I=1,KI 00354 ZTS(I) = PTS(I) 00355 ZTP(I) = PTP(I) 00356 ENDDO 00357 ! 00358 ! 00359 ! Soil analysis based on optimal interpolation 00360 ! 00361 write(*,*) 'PERFORMING OI SOIL ANALYSIS' 00362 CALL OI_CACSTS(KI,ZT2INC,ZH2INC,ZWGINC,PWS_O, & 00363 IDAT,NSSSSS, & 00364 ZTP,ZWP,ZTL,ZSNS,ZTS,ZWS, & 00365 ZTCLS,ZHCLS,ZUCLS,ZVCLS,PSSTC,PWPINC1,PWPINC2,PWPINC3, & 00366 PT2MBIAS,PH2MBIAS, & 00367 PRRCL,PRRSL,PRRCN,PRRSN,PATMNEB,ZEVAP,ZEVAPTR, & 00368 PITM,PVEG,PALBF,PEMISF,PZ0F, & 00369 PIVEG,PARG,PD2,PSAB,PLAI,PRSMIN,PZ0H, & 00370 PTSC,PTPC,PWSC,PWPC,PSNC, & 00371 PGELAT,PGELAM,PGEMU) 00372 00373 ! Update snow 00374 PSWE=ZSNS 00375 00376 ! 00377 ! Perform soil moiture analyses 00378 ! 00379 ZWSINC(:) = 0.0 00380 ZWPINC(:) = 0.0 00381 ZTLINC(:) = 0.0 00382 ! 00383 ZWSINC(:) = ZWS(:) - PWS(:)*(RD1*XRHOLW) 00384 ZWPINC(:) = ZWP(:) - PWP(:)*(PD2(:)*XRHOLW) 00385 ZTLINC(:) = ZTL(:) - PTL(:)*(PD2(:)*XRHOLW) 00386 ! 00387 PWS(:) = ZWS(:)/(RD1*XRHOLW) 00388 PWP(:) = ZWP(:)/(PD2(:)*XRHOLW) 00389 PTL(:) = ZTL(:)/(PD2(:)*XRHOLW) 00390 ! 00391 ! Perform temperature analyses 00392 ! 00393 ! 00394 ZTSINC(:) = 0.0 00395 ZTPINC(:) = 0.0 00396 ! 00397 ZTSINC(:) = ZTS(:) - PTS(:) 00398 ZTPINC(:) = ZTP(:) - PTP(:) 00399 ! 00400 PTS(:) = ZTS(:) 00401 PTP(:) = ZTP(:) 00402 00403 ! 00404 ! 00405 ! PRINT statistics of the soil analysis 00406 ! 00407 PRINT *,'---------------------------------------------------------------' 00408 PRINT *,'Mean WS increments over NATURE ',SUM(ZWSINC)/KI 00409 PRINT *,'Mean WP increments over NATURE ',SUM(ZWPINC)/KI 00410 PRINT *,'Mean TS increments over NATURE ',SUM(ZTSINC)/KI 00411 PRINT *,'Mean TP increments over NATURE ',SUM(ZTPINC)/KI 00412 PRINT *,'Mean TL increments over NATURE ',SUM(ZTLINC)/KI 00413 PRINT *,'---------------------------------------------------------------' 00414 ! 00415 ! Write analysis in LFI file PREP 00416 ! 00417 #ifdef LFI 00418 CFILEOUT_LFI=CPREPFILE 00419 #endif 00420 CALL FLAG_UPDATE(.FALSE.,.TRUE.,.FALSE.,.FALSE.) 00421 CALL INIT_IO_SURF_n(YPROGRAM,'NATURE','SURF ','WRITE') 00422 ! 00423 YVAR='WG1' 00424 YPREFIX='X_Y_WG1 (m3/m3) ' 00425 CALL WRITE_SURF(YPROGRAM,YVAR,PWS,IRESP,HCOMMENT=YPREFIX) 00426 YVAR='WG2' 00427 YPREFIX='X_Y_WG2 (m3/m3) ' 00428 CALL WRITE_SURF(YPROGRAM,YVAR,PWP,IRESP,HCOMMENT=YPREFIX) 00429 YVAR='WGI2' 00430 YPREFIX='X_Y_WGI2 (m3/m3) ' 00431 CALL WRITE_SURF(YPROGRAM,YVAR,PTL,IRESP,HCOMMENT=YPREFIX) 00432 YVAR='TG1' 00433 YPREFIX='X_Y_TG1 (K) ' 00434 CALL WRITE_SURF(YPROGRAM,YVAR,PTS,IRESP,HCOMMENT=YPREFIX) 00435 YVAR='TG2' 00436 YPREFIX='X_Y_TG2 (K) ' 00437 CALL WRITE_SURF(YPROGRAM,YVAR,PTP,IRESP,HCOMMENT=YPREFIX) 00438 ! 00439 CALL END_IO_SURF_n(YPROGRAM) 00440 CALL IO_BUFF_CLEAN_n 00441 00442 ! 00443 ! ------------------------------------------------------------------------------------- 00444 IF (LHOOK) CALL DR_HOOK('ASSIM_NATURE_ISBA_OI',1,ZHOOK_HANDLE) 00445 END SUBROUTINE ASSIM_NATURE_ISBA_OI
1.8.0
