|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 SUBROUTINE ASSIM_NATURE_ISBA_EKF(HPROGRAM, KI, & 00002 PT2M, PHU2M,& 00003 HTEST ) 00004 00005 ! ----------------------------------------------------------------------------- 00006 ! 00007 ! Land Data Assimilation System based on an Extended Kalman Filter 00008 ! 00009 ! Revised version : JFM (15 September 2008) 00010 ! 00011 ! The control vector can be any element of (TG1,TG2,WG1,WG2) - Choice in namelist 00012 ! 00013 ! The observations can be any element of (T2M,HU2M,WG1) - Choice in namelist 00014 ! 00015 ! Possibility to evolve the B matrix in the cycling - otherwise SEKF 00016 ! 00017 ! First version including patches (15 October 2008) 00018 ! 00019 ! ----------------------------------------------------------------------------- 00020 ! 00021 USE MODD_TYPE_DATE_SURF,ONLY : DATE_TIME 00022 USE MODD_ASSIM, ONLY : LPRT,LSIM,LBEV,LBFIXED,NOBSTYPE,YERROBS,INCO,& 00023 IVAR,NVAR,NVARMAX,TPRT_M,XSIGMA_M,XVAR_M,PREFIX_M, & 00024 INCV,SCALE_Q 00025 USE MODN_IO_OFFLINE, ONLY : CPGDFILE,CPREPFILE 00026 ! 00027 #ifdef LFI 00028 USE MODD_IO_SURF_LFI, ONLY : CFILEIN_LFI, CFILEOUT_LFI 00029 #endif 00030 ! 00031 USE MODD_SURF_ATM_n, ONLY : NSIZE_NATURE,NDIM_FULL 00032 ! 00033 USE YOMHOOK, ONLY : LHOOK,DR_HOOK 00034 USE PARKIND1, ONLY : JPRB 00035 ! 00036 USE MODI_ABOR1_SFX 00037 USE MODI_INIT_IO_SURF_n 00038 USE MODI_READ_SURF 00039 USE MODI_END_IO_SURF_n 00040 USE MODI_IO_BUFF_CLEAN_n 00041 USE MODI_WRITE_SURF 00042 USE MODI_ADD_FORECAST_TO_DATE_SURF 00043 ! 00044 ! ----------------------------------------------------------- 00045 ! 00046 IMPLICIT NONE 00047 CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes 00048 INTEGER, INTENT(IN) :: KI 00049 REAL, DIMENSION(KI), INTENT(IN) :: PT2M 00050 REAL, DIMENSION(KI), INTENT(IN) :: PHU2M 00051 CHARACTER(LEN=2), INTENT(IN) :: HTEST ! must be equal to 'OK' 00052 00053 INTEGER :: NOBS 00054 CHARACTER(LEN=28) :: YNAMELIST = 'OPTIONS.nam ' 00055 ! 00056 ! Declarations of local variables 00057 ! 00058 CHARACTER(LEN=3), PARAMETER :: YINIT = 'ALL' 00059 CHARACTER(LEN=6) :: YPROGRAM = 'LFI ' 00060 CHARACTER(LEN=6) :: YPROGRAM2 = 'FA ' 00061 INTEGER :: OBSCOUNT 00062 CHARACTER(LEN=1) :: LCHAR 00063 INTEGER :: IYEAR ! current year (UTC) 00064 INTEGER :: IMONTH ! current month (UTC) 00065 INTEGER :: IDAY ! current day (UTC) 00066 INTEGER :: IHOUR ! current hour (UTC) 00067 REAL :: ZTIME ! current time since start of the run (s) 00068 REAL :: PTSTEP_OUTPUT,PTSTEP_FORC ! OUPUT step, Duration and FORCING Step 00069 INTEGER :: IRESP,PATCH_NUMBER ! return code 00070 INTEGER :: ILOBS ! Namelist unit number 00071 TYPE (DATE_TIME) :: TTIME ! Current date and time 00072 INTEGER :: NBOUTPUT ! Number of time step 00073 INTEGER :: ISTEP ! 00074 REAL,DIMENSION(:,:,:),ALLOCATABLE :: YF ! Vector of model observations (averaged) 00075 REAL,DIMENSION(:,:,:,:),ALLOCATABLE :: YF_PATCH ! vector of model observations (for each pacth) 00076 REAL,DIMENSION(:,:,:,:),ALLOCATABLE :: XF ! Vector of forecast control variables 00077 REAL,DIMENSION(:,:,:),ALLOCATABLE :: XI ! Vector of control variables (at beginning of timestep) 00078 REAL,DIMENSION(:,:,:,:),ALLOCATABLE :: HO ! Jacobian of observation operator 00079 REAL,DIMENSION(:,:,:,:),ALLOCATABLE :: HOWR ! copy of HO for writing out 00080 REAL,DIMENSION(:,:,:,:),ALLOCATABLE :: HOT ! Transpose of HO 00081 REAL,DIMENSION(:,:,:),ALLOCATABLE :: R ! covariance matrix of observation errors 00082 REAL,DIMENSION(:,:,:,:),ALLOCATABLE :: B ! background error covariance matrix 00083 REAL,DIMENSION(:,:),ALLOCATABLE :: IDENT ! identitiy matrix, used for Ba 00084 REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: GAIN ! Kalman gain (used expicitly for Ba) 00085 REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: LTM ! linear tangent matrix for teh f'ward model 00086 REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: Q ! model error matrix 00087 REAL,DIMENSION(:,:),ALLOCATABLE :: YO ! vector of observations 00088 REAL,DIMENSION(:),ALLOCATABLE :: ZCLAY ! Pourcentage of clay (varies from 0 to 1) 00089 REAL,DIMENSION(:),ALLOCATABLE :: ZSAND 00090 REAL,DIMENSION(:),ALLOCATABLE :: COFSWI ! The difference (Wfc - Wwilt) 00091 REAL,DIMENSION(:),ALLOCATABLE :: SMSAT ! 00092 REAL,DIMENSION(:,:,:),ALLOCATABLE :: ZEPS ! The perturbation amplitude 00093 REAL,DIMENSION(:,:,:),ALLOCATABLE :: XINCR ! Analysis increment 00094 REAL,DIMENSION(:,:,:),ALLOCATABLE :: SIMOBS ! Simulated temperature and relative humidity (available per patch ?) 00095 REAL, DIMENSION (:,:),ALLOCATABLE :: XPATCH ! Fraction covered by each patch 00096 REAL,DIMENSION(:,:,:),ALLOCATABLE :: SIMMOD ! Control variables (for B propagation) 00097 REAL,DIMENSION(:,:),ALLOCATABLE :: VECT ! The analysed variable 00098 CHARACTER(LEN=200) :: NMFILE_CANARI ! Name of the observation, perturbed or reference file 00099 CHARACTER(LEN=9) :: HFNAME 00100 CHARACTER(LEN=17) :: LFNAME 00101 INTEGER :: IND 00102 INTEGER :: LTEST 00103 ! 00104 ! Local Matrix for Analysis calculation 00105 ! 00106 REAL,DIMENSION(:,:,:),ALLOCATABLE :: K1 00107 REAL,DIMENSION(:,:),ALLOCATABLE :: ZX,ZB,ZP 00108 REAL,DIMENSION(:,:),ALLOCATABLE :: YOWR 00109 INTEGER :: I,J,K,KK,L 00110 ! 00111 REAL,DIMENSION(:),ALLOCATABLE :: TPRT ! The perturbation amplitude 00112 REAL,DIMENSION(:),ALLOCATABLE :: XSIGMA ! covariance of background errors if B is fixed 00113 ! ! covariance of model errors if B evolving 00114 CHARACTER(LEN=3),DIMENSION(:),ALLOCATABLE :: XVAR ! X is ctrl ! Name of control variables (syntax of surfex in PREP.txt file ) 00115 CHARACTER(LEN=100),DIMENSION(:),ALLOCATABLE:: PREFIX ! The prefix of the control variables (in PREP.txt file) 00116 CHARACTER(LEN=10),DIMENSION(:),ALLOCATABLE :: XOBS ! Identifier for simulated observations 00117 CHARACTER(LEN=3) :: YREAD 00118 ! 00119 INTEGER :: NDIM 00120 00121 INTEGER :: ILUOUT ! ascii output unit number 00122 INTEGER :: ILUNAM ! namelist unit number 00123 INTEGER :: ISTAT 00124 LOGICAL :: GFOUND ! return logical when reading namelist 00125 00126 REAL,DIMENSION(:),ALLOCATABLE :: PT2M_O 00127 REAL,DIMENSION(:),ALLOCATABLE :: PHU2M_O 00128 REAL,DIMENSION(:),ALLOCATABLE :: PEPES 00129 00130 INTEGER :: COMPT 00131 INTEGER :: IVERSION, IBUGFIX 00132 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00133 00134 IF (LHOOK) CALL DR_HOOK('ASSIM_NATURE_ISBA_EKF',0,ZHOOK_HANDLE) 00135 00136 IF (HTEST/='OK') THEN 00137 CALL ABOR1_SFX('ASSIM_NATURE_ISBA_EKF: FATAL ERROR DURING ARGUMENT TRANSFER') 00138 END IF 00139 00140 PRINT * 00141 PRINT *,' --------------------------' 00142 PRINT *,' | ENTERING VARASSIM |' 00143 PRINT *,' --------------------------' 00144 PRINT * 00145 ! 00146 ! 0.2 Allocate arrays depending on the control vector dimension 00147 ! 00148 ALLOCATE (TPRT(NVAR)) 00149 ALLOCATE (XSIGMA(NVAR)) 00150 ALLOCATE (XVAR(NVAR)) 00151 ALLOCATE (PREFIX(NVAR)) 00152 ! 00153 IF (SUM(INCV) /= NVAR) THEN 00154 PRINT *,' INCONSISTENCY in set-up of CONTROL VARIABLES',SUM(INCV),NVAR 00155 STOP 00156 ENDIF 00157 ! 00158 J = 1 00159 DO I = 1,NVARMAX 00160 IF (INCV(I) == 1 .AND. J <= NVAR ) THEN 00161 TPRT(J) = TPRT_M(I) 00162 XSIGMA(J) = XSIGMA_M(I) 00163 XVAR(J) = XVAR_M(I) 00164 PREFIX(J) = PREFIX_M(I) 00165 J = J + 1 00166 ENDIF 00167 ENDDO 00168 ! 00169 ! 0.3 Allocate arrays depending on the number of observation type 00170 ! 00171 ALLOCATE (XOBS(NOBSTYPE)) 00172 ! 00173 ! 1. Initializations 00174 ! 00175 ILOBS=55 00176 ! 00177 ! 00178 ! LFI file handling 00179 ! 00180 #ifdef LFI 00181 CFILEIN_LFI = CPREPFILE ! output of PREP 00182 #endif 00183 ! 00184 ! Find current time 00185 ! 00186 CALL INIT_IO_SURF_n(YPROGRAM,'FULL ','SURF ','READ ') 00187 CALL READ_SURF(YPROGRAM,'DTCUR',TTIME,IRESP) 00188 CALL END_IO_SURF_n(YPROGRAM) 00189 CALL IO_BUFF_CLEAN_n 00190 ! 00191 ! Read number of patches 00192 ! 00193 #ifdef LFI 00194 CFILEIN_LFI = CPGDFILE ! output of PGD 00195 #endif 00196 CALL INIT_IO_SURF_n(YPROGRAM,'NATURE','ISBA ','READ ') 00197 CALL READ_SURF(YPROGRAM,'PATCH_NUMBER',PATCH_NUMBER,IRESP) 00198 CALL READ_SURF(YPROGRAM,'VERSION',IVERSION,IRESP) 00199 CALL READ_SURF(YPROGRAM,'BUG',IBUGFIX,IRESP) 00200 CALL END_IO_SURF_n(YPROGRAM) 00201 CALL IO_BUFF_CLEAN_n 00202 00203 ! 00204 ! Define prefixes for simulated observations 00205 ! 00206 print *,'number of patches =',PATCH_NUMBER 00207 IF (PATCH_NUMBER > 1) THEN 00208 IF (IVERSION<7 .OR. IVERSION==7 .AND. IBUGFIX<3) THEN 00209 XOBS(1) = 'T2M_PATCH' 00210 XOBS(2) = 'HU2M_PATCH' 00211 ELSE 00212 XOBS(1) = 'T2M_P' 00213 XOBS(2) = 'HU2M_P' 00214 ENDIF 00215 ELSE 00216 XOBS(1) = 'T2M_ISBA' 00217 XOBS(2) = 'HU2M_ISBA' 00218 ENDIF 00219 IF (NOBSTYPE > 2) XOBS(3) = 'WG1' 00220 ! 00221 ALLOCATE(XI(NSIZE_NATURE,PATCH_NUMBER,NVAR)) 00222 ALLOCATE(VECT(NSIZE_NATURE,PATCH_NUMBER)) 00223 ALLOCATE(ZSAND(NSIZE_NATURE)) 00224 ALLOCATE(SMSAT(NSIZE_NATURE)) 00225 ! 00226 ! Read in control variables 00227 ! 00228 #ifdef LFI 00229 CFILEIN_LFI = CPREPFILE ! output of PREP 00230 #endif 00231 CALL INIT_IO_SURF_n(YPROGRAM,'NATURE','ISBA ','READ ') 00232 DO L = 1,NVAR 00233 CALL READ_SURF(YPROGRAM,XVAR(L),XI(:,:,L),IRESP) 00234 PRINT *,XVAR(L),' - initial ',XI(1,1,L) 00235 ENDDO 00236 CALL END_IO_SURF_n(YPROGRAM) 00237 CALL IO_BUFF_CLEAN_n 00238 00239 ! 00240 ALLOCATE(COFSWI(NSIZE_NATURE)) 00241 ALLOCATE(ZCLAY(NSIZE_NATURE)) 00242 ALLOCATE(B(NSIZE_NATURE,PATCH_NUMBER,NVAR,NVAR)) 00243 ! 00244 ! Read CLAY fraction to compute the SWI range (Wfc - Wwilt) 00245 ! (XSIGMA is defined in terms of SWI), need to convert to equivalent v/v 00246 ! using same clay fraction in both layers 00247 ! Read SAND fraction to compute the saturation for conversion of ERS SWI 00248 ! 00249 #ifdef LFI 00250 CFILEIN_LFI = CPGDFILE ! output of PGD 00251 #endif 00252 CALL INIT_IO_SURF_n(YPROGRAM,'NATURE','ISBA ','READ ') 00253 CALL READ_SURF(YPROGRAM,'CLAY',ZCLAY(:),IRESP) 00254 CALL READ_SURF(YPROGRAM,'SAND',ZSAND(:),IRESP) 00255 CALL END_IO_SURF_n(YPROGRAM) 00256 CALL IO_BUFF_CLEAN_n 00257 00258 ! 00259 DO I=1,NSIZE_NATURE 00260 COFSWI(I)=0.001*(89.0467*((100.*ZCLAY(I))**0.3496)-37.1342*((100.*ZCLAY(I))**0.5)) 00261 SMSAT(I)=0.001*(-1.08*100*ZSAND(I)+494.305) 00262 ENDDO 00263 ! 00264 ! Frequency of assimilation cycling and data availability 00265 ! 00266 IHOUR = 6 00267 NBOUTPUT = 1 00268 PTSTEP_OUTPUT = 6 00269 ! 00270 ! Time 00271 ! 00272 IYEAR = TTIME%TDATE%YEAR 00273 IMONTH = TTIME%TDATE%MONTH 00274 IDAY = TTIME%TDATE%DAY 00275 ZTIME = TTIME%TIME 00276 00277 ! 00278 ! ==================================================================== 00279 ! 00280 ! ---------------------- 00281 ! VARASSIM OPTION : LPRT 00282 ! ---------------------- 00283 ! 00284 ! ==================================================================== 00285 ! Compute perturbations for control variables 00286 ! (write to CPREPFILE then exit) 00287 ! Element of the control variable chosen by IVAR in namelist 00288 ! 00289 00290 IF ( LPRT ) THEN 00291 ! read in control variable 00292 #ifdef LFI 00293 CFILEOUT_LFI=CPREPFILE 00294 #endif 00295 DO I=1,NSIZE_NATURE 00296 DO J=1,PATCH_NUMBER 00297 VECT(I,J) = XI(I,J,IVAR) + TPRT(IVAR)*XI(I,J,IVAR) 00298 ENDDO 00299 ENDDO 00300 PRINT *,XVAR(IVAR),' - ptrbd', vect(1,1) 00301 CALL INIT_IO_SURF_n(YPROGRAM,'NATURE','ISBA ','WRITE') 00302 CALL WRITE_SURF(YPROGRAM,XVAR(IVAR),VECT,IRESP,PREFIX(IVAR)) 00303 CALL END_IO_SURF_n(YPROGRAM) 00304 CALL IO_BUFF_CLEAN_n 00305 ! 00306 ! Initialisation of B matrix 00307 ! 00308 B(:,:,:,:) = 0.0 00309 DO L=1,NVAR 00310 DO I=1,NSIZE_NATURE 00311 DO J=1,PATCH_NUMBER 00312 IF (XVAR(L) == 'WG2') B(I,J,L,L) = XSIGMA(L)*XSIGMA(L)*COFSWI(I)*COFSWI(I) 00313 IF (XVAR(L) == 'WG1') B(I,J,L,L) = XSIGMA(L)*XSIGMA(L)*COFSWI(I)*COFSWI(I) 00314 IF (XVAR(L) == 'TG2') B(I,J,L,L) = XSIGMA(L)*XSIGMA(L) 00315 IF (XVAR(L) == 'TG1') B(I,J,L,L) = XSIGMA(L)*XSIGMA(L) 00316 ENDDO 00317 ENDDO 00318 ENDDO 00319 ! 00320 PRINT *,'writing out initial B' 00321 OPEN (unit=111,file='BGROUNDin0',status='new',IOSTAT=istat) 00322 IF (istat .NE. 0) THEN 00323 STOP 'BGROUNDin0 already written' 00324 ELSE 00325 WRITE (111,*) ((B(I,J,:,:),J=1,PATCH_NUMBER),I=1,NSIZE_NATURE) 00326 CLOSE(111) 00327 ENDIF 00328 ! 00329 PRINT * 00330 PRINT *,' -----------------------------------' 00331 PRINT *,' | EXITING VARASSIM AFTER LPRT |' 00332 PRINT *,' -----------------------------------' 00333 PRINT * 00334 STOP 00335 ENDIF 00336 00337 00338 ! ==================================================================== 00339 ! 00340 ! ---------------------- 00341 ! VARASSIM OPTION : LSIM 00342 ! ---------------------- 00343 ! 00344 ! ==================================================================== 00345 ! Store simulated observations and evolved prognostic variables. 00346 ! This option is called (NVAR+1) times: 00347 ! 1) with reference initial state (once) 00348 ! 2) with perturbed initial states (NVAR times) 00349 ! 00350 00351 IF ( LSIM ) THEN 00352 ! 00353 #ifdef LFI 00354 CFILEIN_LFI = CPREPFILE ! output of PREP 00355 #endif 00356 CALL INIT_IO_SURF_n(YPROGRAM,'NATURE','ISBA ','READ ') 00357 ! 00358 ALLOCATE(SIMOBS(NSIZE_NATURE,PATCH_NUMBER,NOBSTYPE)) 00359 ALLOCATE(SIMMOD(NSIZE_NATURE,PATCH_NUMBER,NVAR)) 00360 ! 00361 00362 00363 DO I = 1,NOBSTYPE 00364 CALL READ_SURF(YPROGRAM,XOBS(I),SIMOBS(:,:,I),IRESP) 00365 ENDDO 00366 ! 00367 00368 OPEN (unit=111,file='OBSIMU',status='unknown') 00369 DO I=1,NSIZE_NATURE 00370 DO J=1,PATCH_NUMBER 00371 WRITE (111,*) (SIMOBS(I,J,K),K = 1,NOBSTYPE) 00372 ENDDO 00373 ENDDO 00374 CLOSE(111) 00375 ! 00376 DO I = 1, NVAR 00377 CALL READ_SURF(YPROGRAM,XVAR(I),SIMMOD(:,:,I),IRESP) 00378 ENDDO 00379 CALL END_IO_SURF_n(YPROGRAM) 00380 CALL IO_BUFF_CLEAN_n 00381 00382 OPEN (unit=111,file='MDSIMU',status='unknown',IOSTAT=istat) 00383 DO I=1,NSIZE_NATURE 00384 DO J=1,PATCH_NUMBER 00385 WRITE (111,*) (SIMMOD(I,J,L),L=1,NVAR) 00386 ENDDO 00387 ENDDO 00388 CLOSE(111) 00389 ! 00390 PRINT * 00391 PRINT *,' -----------------------------------' 00392 PRINT *,' | EXITING VARASSIM AFTER LSIM |' 00393 PRINT *,' -----------------------------------' 00394 PRINT * 00395 STOP 00396 ENDIF 00397 ! ==================================================================== 00398 ! 00399 ! ---------------------- 00400 ! VARASSIM OPTION : LBEV 00401 ! ---------------------- 00402 ! 00403 ! ==================================================================== 00404 ! Calculate the LTM, and evolve B. 00405 ! 00406 IF (LBEV) THEN 00407 ! 00408 ALLOCATE(LTM(NSIZE_NATURE,PATCH_NUMBER,NVAR,NVAR)) 00409 ALLOCATE(ZEPS(NSIZE_NATURE,PATCH_NUMBER,NVAR)) 00410 ALLOCATE(XF(NSIZE_NATURE,PATCH_NUMBER,NVAR+1,NVAR)) 00411 00412 ! 00413 PRINT *, 'evolving B to time ', IYEAR,IMONTH,IDAY,IHOUR 00414 ! 00415 OPEN (unit=111,file='BGROUNDin',status='unknown',IOSTAT=istat) 00416 READ (111,*) ((B(I,J,:,:),J=1,PATCH_NUMBER),I=1,NSIZE_NATURE) 00417 CLOSE(111) 00418 print *,'read previous B matrix ==>',B(1,1,1,1),NVAR 00419 ! 00420 ! Calculate the TLM of the forecast model 00421 ! 00422 ! a) read in perturbed forecasts 00423 ! 00424 DO L=1,NVAR 00425 WRITE(LCHAR,'(I1) ') L 00426 NMFILE_CANARI='MDSIMU_PERT_'//LCHAR//'_' 00427 CALL GET_FILE_NAME(IYEAR,IMONTH,IDAY,IHOUR,NMFILE_CANARI) 00428 PRINT *,'--> reading in ptb file: ',NMFILE_CANARI 00429 OPEN(UNIT=111,FILE=NMFILE_CANARI,FORM='FORMATTED',STATUS='OLD') 00430 DO I=1,NSIZE_NATURE 00431 DO J=1,PATCH_NUMBER 00432 READ(111,*) (XF(I,J,L+1,K),K=1,NVAR) 00433 ENDDO 00434 ENDDO 00435 CLOSE(111) 00436 PRINT *, 'read in ptbd forecasts for L ',L, XF(1,1,L+1,1) 00437 ENDDO 00438 ! 00439 ! b) read in reference forecasts 00440 ! 00441 NMFILE_CANARI='MDSIMU_REFR_' 00442 PRINT *,'--> reading in reference file' 00443 CALL GET_FILE_NAME(IYEAR,IMONTH,IDAY,IHOUR,NMFILE_CANARI) 00444 OPEN(UNIT=111,FILE=NMFILE_CANARI,FORM='FORMATTED',STATUS='OLD') 00445 ! 00446 DO I=1,NSIZE_NATURE 00447 DO J=1,PATCH_NUMBER 00448 READ(111,*) (XF(I,J,1,K),K=1,NVAR) 00449 ENDDO 00450 ENDDO 00451 CLOSE(111) 00452 PRINT *, 'read in ref forecasts for L',L, XF(1,1,1,1) 00453 ! 00454 ! c) calculate initial perturbation 00455 ! 00456 DO L=1,NVAR 00457 WHERE(XI(:,:,L).NE.999.0) 00458 ZEPS(:,:,L) = TPRT(L)*XI(:,:,L) 00459 ENDWHERE 00460 ENDDO 00461 ! 00462 ! d) calculate LTM - one patch only 00463 ! 00464 DO L=1,NVAR ! control variable (x at previous time step) 00465 DO I=1,NSIZE_NATURE 00466 DO J=1,PATCH_NUMBER 00467 DO K=1,NVAR 00468 LTM(I,J,K,L) = (XF(I,J,L+1,K) - XF(I,J,1,K))/ZEPS(I,J,L) ! Jacobian of fwd model 00469 ENDDO 00470 ENDDO 00471 ENDDO 00472 ENDDO 00473 ! 00474 ! evolve B - one patch only 00475 ! 00476 DO I=1,NSIZE_NATURE 00477 DO J=1,PATCH_NUMBER 00478 B(I,J,:,:) = MATMUL(LTM(I,J,:,:),MATMUL(B(I,J,:,:),TRANSPOSE(LTM(I,J,:,:)))) 00479 ENDDO 00480 ENDDO 00481 ! 00482 PRINT *,'LTM d(wg2)/d(wg2)', LTM(1,1,1,1) 00483 ! 00484 ! write out the LTM for the forward model 00485 ! 00486 DO L=1,NVAR 00487 DO K=1,NVAR 00488 WRITE(LCHAR,'(I1)') K 00489 LFNAME='LTM_del'//XVAR(K)//'_del'//XVAR(L) 00490 OPEN(UNIT=111,FILE=LFNAME,FORM='FORMATTED',STATUS='UNKNOWN',POSITION='APPEND') 00491 DO I=1,NSIZE_NATURE 00492 DO J=1,PATCH_NUMBER 00493 WRITE (111,*) LTM(I,J,K,L) 00494 ENDDO 00495 ENDDO 00496 CLOSE(111) 00497 ENDDO 00498 ENDDO 00499 ! 00500 ! write out current B (for use in next cycle) 00501 ! 00502 print *,'store B matrix after TL evolution ==>',B(1,1,1,1) 00503 PRINT *,'writing out B' 00504 OPEN (unit=111,file='BGROUNDout',status='unknown') 00505 WRITE (111,*) ((B(I,J,:,:),J=1,PATCH_NUMBER),I=1,NSIZE_NATURE) 00506 CLOSE(111) 00507 ! 00508 DEALLOCATE(B) 00509 DEALLOCATE(LTM) 00510 DEALLOCATE(ZEPS) 00511 DEALLOCATE(XF) 00512 ! 00513 PRINT * 00514 PRINT *,' -----------------------------------' 00515 PRINT *,' | EXITING VARASSIM AFTER LBEV |' 00516 PRINT *,' -----------------------------------' 00517 PRINT * 00518 STOP 00519 ENDIF 00520 ! ==================================================================== 00521 ! 00522 ! if not LSIM,LPTR, or LBEV proceed with analysis 00523 ! 00524 ! Number of available observation files 00525 PRINT *,'--> through to analysis section' 00526 00527 00528 00529 ! NOBS=2 00530 00531 00532 ! 00533 ! NOBS=0 00534 !R ZTIME=PTSTEP_OUTPUT 00535 ! 00536 ! DO ISTEP=1,NBOUTPUT 00537 !R CALL ADD_FORECAST_TO_DATE_SURF(IYEAR,IMONTH,IDAY,ZTIME) 00538 !R ZTIME = ZTIME + PTSTEP_OUTPUT 00539 !R NMFILE_CANARI='CANARI_NATURE_' 00540 !R CALL GET_FILE_NAME(IYEAR,IMONTH,IDAY,IHOUR,NMFILE_CANARI) 00541 !R print *,'observation file name =>',NMFILE_CANARI 00542 !R OPEN(UNIT=ILOBS,FILE=NMFILE_CANARI,FORM='UNFORMATTED',STATUS='OLD',ERR=22) 00543 ! NOBS=NOBS+NOBSTYPE 00544 ! 22 CONTINUE 00545 !R CLOSE(ILOBS) 00546 ! ENDDO 00547 ! 00548 !R IF (NOBS == 0) THEN 00549 !R PRINT *,'No observations read in OBS file - stop' 00550 !R STOP 00551 !R ENDIF 00552 ! 00553 00554 00555 00556 00557 ALLOCATE (PT2M_O(NSIZE_NATURE)) 00558 ALLOCATE (PHU2M_O(NSIZE_NATURE)) 00559 00560 COMPT=0 00561 DO I=1,KI 00562 COMPT=COMPT+1 00563 PT2M_O(COMPT)=PT2M(I) 00564 PHU2M_O(COMPT)=PHU2M(I) 00565 ENDDO 00566 00567 00568 ! 00569 ! Time reinitialization 00570 ! 00571 IYEAR = TTIME%TDATE%YEAR 00572 IMONTH = TTIME%TDATE%MONTH 00573 IDAY = TTIME%TDATE%DAY 00574 ZTIME = TTIME%TIME 00575 ! 00576 ! Allocation 00577 ! 00578 ! Perturbed simulations 00579 ! 00580 ALLOCATE(YF(NSIZE_NATURE,NVAR+1,NOBSTYPE)) 00581 ALLOCATE(YF_PATCH(NSIZE_NATURE,PATCH_NUMBER,NVAR+1,NOBSTYPE)) 00582 ALLOCATE(ZEPS(NSIZE_NATURE,PATCH_NUMBER,NVAR)) 00583 ALLOCATE(XF(NSIZE_NATURE,PATCH_NUMBER,NVAR+1,NVAR)) 00584 ALLOCATE(XPATCH(NSIZE_NATURE,PATCH_NUMBER)) 00585 ! 00586 ! Read fraction of each patch (=> LPGD=T in OPTIONS.nam) 00587 ! 00588 #ifdef LFI 00589 CFILEIN_LFI = CPGDFILE ! output of PGD 00590 #endif 00591 CALL INIT_IO_SURF_n(YPROGRAM,'NATURE','ISBA ','READ ') 00592 IF (PATCH_NUMBER > 1) THEN 00593 CALL READ_SURF(YPROGRAM,'PATCH',XPATCH,IRESP) 00594 ELSE 00595 XPATCH(:,1) = 1.0 00596 ENDIF 00597 CALL END_IO_SURF_n(YPROGRAM) 00598 CALL IO_BUFF_CLEAN_n 00599 00600 ! 00601 00602 00603 ! Initial values (to be analysed) 00604 ! Observations 00605 ! 00606 ALLOCATE(YO(NSIZE_NATURE,NOBSTYPE)) 00607 ALLOCATE(YOWR(NDIM_FULL,NOBSTYPE)) 00608 ! 00609 ! Temporary vectors used by the EKF approach 00610 ! 00611 ALLOCATE(ZB(NSIZE_NATURE,NOBSTYPE)) 00612 ALLOCATE(HO(NSIZE_NATURE,PATCH_NUMBER,NOBSTYPE,NVAR)) 00613 ALLOCATE(HOWR(NSIZE_NATURE,PATCH_NUMBER,NOBSTYPE,NVAR)) 00614 ALLOCATE(HOT(NSIZE_NATURE,PATCH_NUMBER,NVAR,NOBSTYPE)) 00615 ALLOCATE(R(NSIZE_NATURE,NOBSTYPE,NOBSTYPE)) 00616 ALLOCATE(IDENT(NVAR,NVAR)) 00617 ALLOCATE(GAIN(NSIZE_NATURE,PATCH_NUMBER,NVAR,NOBSTYPE)) 00618 ALLOCATE(ZX(NSIZE_NATURE,NOBSTYPE)) 00619 ALLOCATE(ZP(NSIZE_NATURE,NOBSTYPE)) 00620 ALLOCATE(XINCR(NSIZE_NATURE,PATCH_NUMBER,NVAR)) 00621 ALLOCATE(K1(NSIZE_NATURE,NOBSTYPE,NOBSTYPE)) 00622 ! 00623 ! Initialisations 00624 ! 00625 HO(:,:,:,:) = 999.0 ! Linearized observation matrix 00626 HOWR(:,:,:,:) = 999.0 00627 R(:,:,:) = 0.0 ! Observation error matrix 00628 B(:,:,:,:) = 0.0 ! Background error matrix 00629 YOWR(:,:) = 999.0 ! Observation vector on the full grid to be written on file 00630 ! YO(:,:) = 999.0 00631 ZB(:,:) = 999.0 ! Innovation vector 00632 YF(:,:,:) = 0.0 ! Tile averaged simulated observation vector 00633 OBSCOUNT = 0 00634 ! 00635 IDENT=0 ! identity matrix 00636 DO I = 1,NVAR 00637 IDENT(I,I) = 1 00638 ENDDO 00639 ! 00640 ! Calculate delta for control variables 00641 ! 00642 DO L=1,NVAR 00643 ZEPS(:,:,L) = TPRT(L)*XI(:,:,L) 00644 ENDDO 00645 ! 00646 ! NOBS=2 00647 ZTIME=PTSTEP_OUTPUT 00648 ! 00649 ! BEGINNING OF TIME LOOP 00650 ! 00651 00652 00653 TIMELOOP : DO ISTEP=1,NBOUTPUT 00654 ! 00655 ! Update date 00656 ! 00657 00658 00659 CALL ADD_FORECAST_TO_DATE_SURF(IYEAR, IMONTH, IDAY, ZTIME) 00660 ZTIME = ZTIME + PTSTEP_OUTPUT 00661 !R NMFILE_CANARI='CANARI_NATURE_' 00662 !R CALL GET_FILE_NAME(IYEAR,IMONTH,IDAY,IHOUR,NMFILE_CANARI) 00663 ! 00664 ! Open observation files 00665 ! 00666 !R OPEN(UNIT=ILOBS,FILE=NMFILE_CANARI,FORM='FORMATTED',STATUS='OLD',ERR=10) 00667 ! 00668 ! NOBS=NOBS+NOBSTYPE 00669 ! 00670 ! If it exists, read T2m and HU2m observations 00671 ! 00672 !R DO I=1,NSIZE_NATURE 00673 !R READ (ILOBS,*) (YO(I,J),J=1,NOBSTYPE) 00674 !R IF (YO(I,3) == 110.0) YO(I,3) = 999.0 ! frozen soils 00675 !R ENDDO 00676 00677 00678 00679 YO(:,1)=PT2M_O(:) 00680 YO(:,2)=PHU2M_O(:) 00681 00682 PRINT *,'read in obs: ', YO(1,1), YO(1,2), NOBS 00683 00684 00685 00686 ! 00687 ! Calculate perturbations at the date of observations and innovation 00688 ! (2nd dim of YF: index 1 is reference, 2,3... perturbed by cntrl var 1,2... 00689 ! 00690 DO L=1,NVAR 00691 WRITE(LCHAR,'(I1) ') L 00692 NMFILE_CANARI='OBSIMU_PERT_'//LCHAR//'_' 00693 CALL GET_FILE_NAME(IYEAR,IMONTH,IDAY,IHOUR,NMFILE_CANARI) 00694 PRINT *,'--> reading in ptb file: ',NMFILE_CANARI 00695 OPEN(UNIT=111,FILE=NMFILE_CANARI,FORM='FORMATTED',STATUS='OLD') 00696 DO I=1,NSIZE_NATURE 00697 DO J=1,PATCH_NUMBER 00698 ! READ(111,*) (YF_PATCH(I,J,L+1,K+NOBS-NOBSTYPE),K=1,NOBSTYPE) 00699 READ(111,*) (YF_PATCH(I,J,L+1,K),K=1,NOBSTYPE) 00700 ENDDO 00701 ENDDO 00702 CLOSE(111) 00703 PRINT *, 'read in ptbd forecasts for H',L, YF_PATCH(1,1,L+1,1), YF_PATCH(1,1,L+1,2) 00704 ENDDO 00705 00706 00707 00708 ! 00709 NMFILE_CANARI='OBSIMU_REFR_' 00710 PRINT *,'--> reading in reference file' 00711 CALL GET_FILE_NAME(IYEAR,IMONTH,IDAY,IHOUR,NMFILE_CANARI) 00712 OPEN(UNIT=111,FILE=NMFILE_CANARI,FORM='FORMATTED',STATUS='OLD') 00713 ! 00714 DO I=1,NSIZE_NATURE 00715 DO J=1,PATCH_NUMBER 00716 ! READ(111,*) (YF_PATCH(I,J,1,K+NOBS-NOBSTYPE),K=1,NOBSTYPE) 00717 READ(111,*) (YF_PATCH(I,J,1,K),K=1,NOBSTYPE) 00718 ENDDO 00719 ENDDO 00720 CLOSE(111) 00721 PRINT *, 'read in ref forecasts for', YF_PATCH(1,1,1,1),YF_PATCH(1,1,1,2) 00722 10 CONTINUE ! if T2m HU2m observations does not exist 00723 CLOSE(ILOBS) 00724 00725 00726 00727 ! 00728 ! Mean simulated obs averaged over tiles 00729 ! 00730 DO I=1,NSIZE_NATURE 00731 DO J=1,PATCH_NUMBER 00732 IF (XPATCH(I,J) > 0.0) THEN 00733 YF(I,:,:) = YF(I,:,:) + XPATCH(I,J)*YF_PATCH(I,J,:,:) 00734 ENDIF 00735 ENDDO 00736 ENDDO 00737 ! 00738 ENDDO TIMELOOP 00739 00740 00741 00742 ! 00743 ! Rescale the satellite derived soil moisture observations 00744 ! 00745 !R DO I=1,NSIZE_NATURE 00746 !R IF (YO(I,3) /= 999.0) THEN 00747 ! CDF matching on ERS (revised formula) 00748 ! YO(I,3) = (1.16117E-7*YO(I,3)**5 - 2.85771E-5*YO(I,3)**4 + 0.00235291*YO(I,3)**3 & 00749 ! & - 0.070169*YO(I,3)**2 + 1.261382*YO(I,3) + 12.2607)*SMSAT(I)/100. 00750 ! CDF matching for ASCAT soil moisture 00751 !R YO(I,3) = (8.80461E-8*YO(I,3)**5 - 2.21598E-5*YO(I,3)**4 + 0.00188043*YO(I,3)**3 & 00752 !R & - 0.0575883*YO(I,3)**2 + 1.0249301*YO(I,3) + 15.7502)*SMSAT(I)/100. 00753 !R ENDIF 00754 !R ENDDO 00755 ! 00756 00757 00758 00759 00760 ! SET OBSERVATION ERROR 00761 ! 00762 DO I=1,NSIZE_NATURE 00763 DO K=1,NOBSTYPE 00764 IF (K .LT. 3) THEN 00765 R(I,K,K) = YERROBS(K)*YERROBS(K) 00766 ELSE 00767 ! convert R for wg1 from SWI to abs value 00768 R(I,K,K) = YERROBS(K)*YERROBS(K)*COFSWI(I)*COFSWI(I) 00769 ENDIF 00770 ENDDO 00771 ENDDO 00772 ! 00773 ! WRITE OUT OBS AND YERROR FOR DIAGNOSTIC PURPOSES 00774 ! 00775 OPEN (unit=111,file='OBSERRORout',status='unknown',IOSTAT=istat) 00776 WRITE (111,*) (R(I,2,2), I=1,NSIZE_NATURE) 00777 CLOSE(111) 00778 OPEN (unit=111,file='OBSout',status='unknown',IOSTAT=istat) 00779 WRITE (111,*) (YO(I,2), I=1,NSIZE_NATURE) 00780 CLOSE(111) 00781 00782 00783 00784 00785 !-------------------------------------------------------------------- 00786 ! 00787 ! Background error matrix 00788 ! 00789 !-------------------------------------------------------------------- 00790 IF (LBFIXED) THEN 00791 DO L=1,NVAR 00792 DO I=1,NSIZE_NATURE 00793 DO J=1,PATCH_NUMBER 00794 IF (XVAR(L) == 'WG2') B(I,J,L,L) = XSIGMA(L)*XSIGMA(L)*COFSWI(I)*COFSWI(I) 00795 IF (XVAR(L) == 'WG1') B(I,J,L,L) = XSIGMA(L)*XSIGMA(L)*COFSWI(I)*COFSWI(I) 00796 IF (XVAR(L) == 'TG2') B(I,J,L,L) = XSIGMA(L)*XSIGMA(L) 00797 IF (XVAR(L) == 'TG1') B(I,J,L,L) = XSIGMA(L)*XSIGMA(L) 00798 ENDDO 00799 ENDDO 00800 ENDDO 00801 ELSE 00802 ALLOCATE(Q(NSIZE_NATURE,PATCH_NUMBER,NVAR,NVAR)) 00803 Q(:,:,:,:) = 0.0 00804 DO L=1,NVAR 00805 DO I=1,NSIZE_NATURE 00806 DO J=1,PATCH_NUMBER 00807 IF (XVAR(L) == 'WG2') Q(I,J,L,L) = SCALE_Q*SCALE_Q*XSIGMA(L)*XSIGMA(L)*COFSWI(I)*COFSWI(I) 00808 IF (XVAR(L) == 'WG1') Q(I,J,L,L) = SCALE_Q*SCALE_Q*XSIGMA(L)*XSIGMA(L)*COFSWI(I)*COFSWI(I) 00809 IF (XVAR(L) == 'TG2') Q(I,J,L,L) = SCALE_Q*SCALE_Q*XSIGMA(L)*XSIGMA(L) 00810 IF (XVAR(L) == 'TG1') Q(I,J,L,L) = SCALE_Q*SCALE_Q*XSIGMA(L)*XSIGMA(L) 00811 ENDDO 00812 ENDDO 00813 ENDDO 00814 PRINT *,'reading in B, then updating with Q' 00815 OPEN (unit=111,file='BGROUNDin',status='unknown',IOSTAT=istat) 00816 READ (111,*) ((B(I,J,:,:),J=1,PATCH_NUMBER),I=1,NSIZE_NATURE) 00817 CLOSE(111) 00818 ! 00819 ! B is the forecast matrix - need to add Q 00820 ! 00821 print *,'B before wg2 wg2 ==> ',sqrt(B(1,1,1,1))/COFSWI(1),B(1,1,1,1) 00822 print *,'Q value wg2 wg2 ==> ',sqrt(Q(1,1,1,1))/COFSWI(1),Q(1,1,1,1) 00823 B = B + Q 00824 print *,'B after wg2 wg2 ==>',sqrt(B(1,1,1,1))/COFSWI(1),B(1,1,1,1) 00825 ! 00826 DEALLOCATE(Q) 00827 ENDIF 00828 ! 00829 00830 00831 00832 ! Data type selection before assimilation (only if NOBS = NOBSTYPE) 00833 ! 00834 IF (NOBSTYPE == NOBSTYPE) THEN 00835 DO I=1,NOBSTYPE 00836 IF (INCO(I) == 0) THEN 00837 YO (:,I) = 999.0 00838 PRINT *,'OBSERVATION TYPE ',XOBS(I),' REMOVED' 00839 ENDIF 00840 ENDDO 00841 ENDIF 00842 ! YO(:,:) = 999.0 ! WARNING this instruction remove all observations 00843 ! 00844 PRINT *,'calculating jacobians',NOBS 00845 DO L=1,NVAR 00846 DO I=1,NSIZE_NATURE 00847 DO J=1,PATCH_NUMBER 00848 DO K=1,NOBSTYPE 00849 HOWR(I,J,K,L) = XPATCH(I,J)*(YF_PATCH(I,J,L+1,K) - YF_PATCH(I,J,1,K))/ZEPS(I,J,L) 00850 IF(YO(I,K) .NE. 999.0) THEN 00851 HO(I,J,K,L) = XPATCH(I,J)*(YF_PATCH(I,J,L+1,K) - YF_PATCH(I,J,1,K))/ZEPS(I,J,L) ! Jacobian of obs operator 00852 ZB(I,K) = YO(I,K) - YF(I,1,K) ! innovation vector 00853 OBSCOUNT = OBSCOUNT + 1 00854 ELSE 00855 ! set obs op and obs innovation to zero if no obs are present 00856 HO(I,J,K,:) = 0.0 00857 ZB(I,K) = 0.0 00858 ENDIF 00859 ENDDO 00860 ENDDO 00861 ENDDO 00862 ENDDO 00863 00864 00865 00866 00867 !----------------------------------------------------- 00868 ! 00869 ! ****** SOIL ANALYSIS ******* 00870 ! 00871 !----------------------------------------------------- 00872 PRINT *,'PERFORMING ANALYSIS' 00873 DO I=1,NSIZE_NATURE 00874 DO J=1,PATCH_NUMBER 00875 HOT(I,J,:,:) = TRANSPOSE(HO(I,J,:,:)) 00876 K1(I,:,:) = MATMUL(HO(I,J,:,:),MATMUL(B(I,J,:,:),HOT(I,J,:,:))) + R(I,:,:) 00877 CALL CHOLDC(NOBSTYPE,K1(I,:,:),ZP(I,:)) ! Cholesky decomposition (1) 00878 CALL CHOLSL(NOBSTYPE,K1(I,:,:),ZP(I,:),ZB(I,:),ZX(I,:)) ! Cholesky decomposition (2) 00879 XINCR(I,J,:) = MATMUL(B(I,J,:,:),MATMUL(HOT(I,J,:,:),ZX(I,:))) 00880 !TA: A possible safety check for WG2 and WG1. Is something needed? I have seen negative values for soil water 00881 DO L=1,NVAR 00882 !IF ( L <= 2 ) THEN 00883 ! XI(I,J,L) = MAX(0.000001,(XI(I,J,L) + XINCR(I,J,L))) 00884 !ELSE 00885 XI(I,J,L) = XI(I,J,L) + XINCR(I,J,L) 00886 !ENDIF 00887 ENDDO 00888 ENDDO 00889 ENDDO 00890 ! 00891 ! *** Write analysis results in ASCII file *** 00892 ! 00893 OPEN (unit=111,file='ANAL_INCR',status='unknown',IOSTAT=istat) 00894 DO I=1,NSIZE_NATURE 00895 WRITE(111,*) (XI(I,1,J), XINCR(I,1,J),J=1,NVAR),ZB(I,1),ZB(I,2) 00896 ENDDO 00897 CLOSE(UNIT=111) 00898 00899 ! ********************************************** 00900 DO L=1,NVAR 00901 PRINT *,'Sum and mean increments for ',TRIM(XVAR(l)),'=',SUM(XINCR(:,1,L)),SUM(XINCR(:,1,L))/REAL(NSIZE_NATURE) 00902 ENDDO 00903 ! 00904 ! Write analysis in LFI file (PREP.lfi) 00905 ! 00906 #ifdef LFI 00907 CFILEOUT_LFI=CPREPFILE 00908 #endif 00909 CALL INIT_IO_SURF_n(YPROGRAM,'NATURE','ISBA ','WRITE') 00910 DO L = 1,NVAR 00911 CALL WRITE_SURF(YPROGRAM,XVAR(L),XI(:,:,L),IRESP,PREFIX(L)) 00912 ENDDO 00913 CALL END_IO_SURF_n(YPROGRAM) 00914 CALL IO_BUFF_CLEAN_n 00915 00916 ! 00917 ! Analysis of B (for use in next cycle) 00918 ! Ba = (I-KH)Bf 00919 ! K = BfHT{K1}**-1 00920 ! K1 needs PATCH dim added 00921 ! 00922 DO I=1,NSIZE_NATURE 00923 DO J=1,PATCH_NUMBER 00924 ! K1 = (R+H.B.HT) (calculate inverse -> output goes to K1) 00925 CALL INVERSE_MATRIX(NOBSTYPE,K1(I,:,:),ZP(I,:)) 00926 GAIN(I,J,:,:) = MATMUL(B(I,J,:,:),MATMUL(HOT(I,J,:,:),K1(I,:,:))) 00927 IF (.NOT.LBFIXED) B(I,J,:,:) = MATMUL((IDENT - MATMUL(GAIN(I,J,:,:),HO(I,J,:,:))),B(I,J,:,:)) 00928 ENDDO 00929 ENDDO 00930 ! 00931 ! Write out analysed B (for use in next cycle) 00932 ! 00933 OPEN (unit=111,file='BGROUNDout',status='unknown',IOSTAT=istat) 00934 WRITE (111,*) ((B(I,J,:,:),J=1,PATCH_NUMBER),I=1,NSIZE_NATURE) 00935 CLOSE(111) 00936 ! 00937 ! **** Write out the observation operator + Gain matrix **** 00938 ! 00939 DO L=1,NVAR 00940 DO K=1, NOBSTYPE 00941 WRITE(LCHAR,'(I1)') K 00942 HFNAME='HO_'//XVAR(L)//'_v'//LCHAR 00943 OPEN(UNIT=111,FILE=HFNAME,FORM='FORMATTED',STATUS='NEW',IOSTAT=istat) 00944 DO I=1,NSIZE_NATURE 00945 DO J=1,PATCH_NUMBER 00946 WRITE(111,*) HOWR(I,J,K,L),GAIN(I,J,L,K) 00947 ENDDO 00948 ENDDO 00949 CLOSE(111) 00950 ENDDO 00951 ENDDO 00952 00953 DEALLOCATE(YOWR) 00954 DEALLOCATE(VECT) 00955 DEALLOCATE(YF) 00956 DEALLOCATE(YF_PATCH) 00957 DEALLOCATE(ZB) 00958 DEALLOCATE(HO) 00959 DEALLOCATE(R) 00960 DEALLOCATE(B) 00961 DEALLOCATE(GAIN) 00962 DEALLOCATE(ZP) 00963 DEALLOCATE(ZEPS) 00964 DEALLOCATE(YO) 00965 DEALLOCATE(HOT) 00966 DEALLOCATE(K1) 00967 DEALLOCATE(ZX) 00968 DEALLOCATE(XINCR) 00969 DEALLOCATE(HOWR) 00970 DEALLOCATE(PT2M_O) 00971 DEALLOCATE(PHU2M_O) 00972 00973 00974 00975 ! 00976 PRINT * 00977 PRINT *,' ---------------------------------------' 00978 PRINT *,' | EXITING VARASSIM AFTER ANALYSIS |' 00979 PRINT *,' ---------------------------------------' 00980 PRINT * 00981 PRINT *,'Number of assimilated observations =',OBSCOUNT 00982 PRINT * 00983 ! 00984 00985 IF (LHOOK) CALL DR_HOOK('ASSIM_NATURE_ISBA_EKF',1,ZHOOK_HANDLE) 00986 00987 CONTAINS 00988 00989 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Subroutine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00990 00991 SUBROUTINE GET_FILE_NAME(IYEAR,IMONTH,IDAY,IHOUR,NMFILE) 00992 USE MODI_TRANS_CHAINE 00993 IMPLICIT NONE 00994 INTEGER, INTENT (IN) :: IYEAR, IMONTH, IDAY, IHOUR 00995 CHARACTER(LEN=50) :: CMONTH,CDAY,CYEAR,CHOUR 00996 CHARACTER(LEN=200), INTENT(INOUT) :: NMFILE 00997 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00998 ! 00999 IF (LHOOK) CALL DR_HOOK('ASSIM_NATURE_ISBA_EKF:GET_FILE_NAME',0,ZHOOK_HANDLE) 01000 ! 01001 CALL TRANS_CHAINE(CMONTH,IMONTH,2) 01002 CALL TRANS_CHAINE(CDAY,IDAY,2) 01003 CALL TRANS_CHAINE(CYEAR,IYEAR,4) 01004 CALL TRANS_CHAINE(CHOUR,IHOUR,2) 01005 IF (IHOUR.eq.0.or.IHOUR.eq.24.or.IHOUR.eq.48) CHOUR='00' 01006 ! 01007 NMFILE = TRIM(NMFILE)//CYEAR(3:4) 01008 NMFILE = TRIM(NMFILE)//CMONTH 01009 NMFILE = TRIM(NMFILE)//CDAY 01010 NMFILE = TRIM(NMFILE)//'H' 01011 NMFILE = TRIM(NMFILE)//CHOUR 01012 NMFILE = TRIM(NMFILE)//'.DAT' 01013 ! 01014 IF (LHOOK) CALL DR_HOOK('ASSIM_NATURE_ISBA_EKF:GET_FILE_NAME',1,ZHOOK_HANDLE) 01015 ! 01016 END SUBROUTINE GET_FILE_NAME 01017 01018 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01019 01020 SUBROUTINE INVERSE_MATRIX(N,A,P) 01021 !-------------------------------------------------------- 01022 ! 01023 ! Explicit inversed matrix after Cholesky decomposition 01024 ! 01025 !-------------------------------------------------------- 01026 IMPLICIT NONE 01027 INTEGER, INTENT(IN) :: N 01028 REAL, DIMENSION (N,N), INTENT(INOUT) :: A 01029 REAL, DIMENSION (N), INTENT(IN) :: P 01030 REAL ZSUM 01031 INTEGER :: I, J, K 01032 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01033 ! 01034 IF (LHOOK) CALL DR_HOOK('ASSIM_NATURE_ISBA_EKF:INVERSE_MATRIX',0,ZHOOK_HANDLE) 01035 ! 01036 DO I=1,N 01037 A(I,I)=1./P(I) 01038 DO J=I+1,N 01039 ZSUM = 0. 01040 DO K=I,J-1 01041 ZSUM = ZSUM - A(J,K)*A(K,I) 01042 ENDDO 01043 A(J,I) = ZSUM/P(J) 01044 ENDDO 01045 ENDDO 01046 DO I=1,N 01047 DO J=I+1,N 01048 A(I,J) =0.0 01049 ENDDO 01050 ENDDO 01051 A = MATMUL(TRANSPOSE(A),A) 01052 ! 01053 IF (LHOOK) CALL DR_HOOK('ASSIM_NATURE_ISBA_EKF:INVERSE_MATRIX',1,ZHOOK_HANDLE) 01054 ! 01055 END SUBROUTINE INVERSE_MATRIX 01056 01057 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01058 01059 SUBROUTINE CHOLDC(N,A,P) 01060 IMPLICIT NONE 01061 INTEGER, INTENT(IN) :: N 01062 REAL, DIMENSION(N,N), INTENT(INOUT) :: A 01063 REAL, DIMENSION(N), INTENT(OUT) :: P 01064 INTEGER :: I 01065 REAL :: SUMM 01066 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01067 ! 01068 IF (LHOOK) CALL DR_HOOK('ASSIM_NATURE_ISBA_EKF:CHOLDC',0,ZHOOK_HANDLE) 01069 ! 01070 DO I=1,N 01071 SUMM = A(I,I)- DOT_PRODUCT(A(I,1:I-1),A(I,1:I-1)) 01072 IF (SUMM <= 0.0) THEN 01073 PRINT*,'ERROR IN CHOLDC' 01074 STOP 01075 ENDIF 01076 P(I) = SQRT(SUMM) 01077 A(I+1:N,I) = (A(I,I+1:N) - MATMUL(A(I+1:N,1:I-1),A(I,1:I-1)))/P(I) 01078 ENDDO 01079 ! 01080 IF (LHOOK) CALL DR_HOOK('ASSIM_NATURE_ISBA_EKF:CHOLDC',1,ZHOOK_HANDLE) 01081 ! 01082 END SUBROUTINE CHOLDC 01083 01084 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01085 01086 SUBROUTINE CHOLSL(N,A,P,B,X) 01087 IMPLICIT NONE 01088 INTEGER, INTENT(IN) :: N 01089 REAL, DIMENSION (N,N), INTENT(IN) :: A 01090 REAL, DIMENSION (N), INTENT(IN) :: P,B 01091 REAL, DIMENSION (N), INTENT(INOUT) :: X 01092 INTEGER :: I 01093 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01094 ! 01095 IF (LHOOK) CALL DR_HOOK('ASSIM_NATURE_ISBA_EKF:CHOLSL',0,ZHOOK_HANDLE) 01096 ! 01097 DO I=1,N 01098 X(I) = (B(I) - DOT_PRODUCT(A(I,1:I-1),X(1:I-1)))/P(I) 01099 ENDDO 01100 DO I=N,1,-1 01101 X(I) = (X(I) - DOT_PRODUCT(A(I+1:N,I),X(I+1:N)))/P(I) 01102 ENDDO 01103 ! 01104 IF (LHOOK) CALL DR_HOOK('ASSIM_NATURE_ISBA_EKF:CHOLSL',1,ZHOOK_HANDLE) 01105 ! 01106 END SUBROUTINE CHOLSL 01107 01108 END SUBROUTINE ASSIM_NATURE_ISBA_EKF 01109
1.8.0
