SURFEX v7.3: /home/dasprezs/EXPORT_v7_3/src/LIB/TOPD/coupl

Go to the documentation of this file.
00001 !-----------------------------------------------------------------
00002 !     #####################
00003       SUBROUTINE COUPL_TOPD(HPROGRAM,HSTEP,KI,KSTEP)
00004 !     #####################
00005 !
00006 !!****  *COUPL_TOPD*  
00007 !!
00008 !!    PURPOSE
00009 !!    -------
00010 !
00011 !     
00012 !         
00013 !     
00014 !!**  METHOD
00015 !!    ------
00016 !
00017 !!    EXTERNAL
00018 !!    --------
00019 !!
00020 !!    none
00021 !!
00022 !!    IMPLICIT ARGUMENTS
00023 !!    ------------------ 
00024 !!
00025 !!    
00026 !!    
00027 !!
00028 !!      
00029 !!    REFERENCE
00030 !!    ---------
00031 !!
00032 !!    
00033 !!      
00034 !!    AUTHOR
00035 !!    ------
00036 !!
00037 !!      K. Chancibault  * LTHE / Meteo-France *
00038 !!
00039 !!    MODIFICATIONS
00040 !!    -------------
00041 !!
00042 !!      Original   15/10/2003
00043 !!      09/2007 : New organisation of exfiltration, computation of saturated
00044 !!                area, routing.
00045 !!                Soil ice content taken into account
00046 !-------------------------------------------------------------------------------
00047 !
00048 !*       0.     DECLARATIONS
00049 !               ------------
00050 !
00051 USE MODD_TOPODYN,        ONLY : NNCAT, NMESHT, NNMC, XMPARA
00052 USE MODD_COUPLING_TOPD,  ONLY : XWG_FULL, XDTOPI, XKAC_PRE, XDTOPT, XWTOPT, XWSTOPT, XAS_NATURE,&
00053                                   XKA_PRE, NMASKT, XWSUPSAT,&
00054                                   XRUNOFF_TOP, XATOP, XWFCTOPI, NNPIX,&
00055                                   XFRAC_D2, XFRAC_D3, XWSTOPI, XDMAXFC, XWFCTOPT, XWGI_FULL,&
00056                                   NFREQ_MAPS_ASAT, XAVG_RUNOFFCM, XAVG_DRAINCM
00057 !
00058 USE MODD_ISBA_n,           ONLY : XWG, XDG, XWGI
00059 USE MODD_CSTS,             ONLY : XRHOLW, XRHOLI
00060 USE MODD_DIAG_EVAP_ISBA_n, ONLY : XAVG_RUNOFFC, XAVG_DRAINC
00061 USE MODD_SURF_ATM_n,       ONLY : NR_NATURE
00062 USE MODD_SURF_ATM_GRID_n,  ONLY : XMESH_SIZE
00063 USE MODD_SURF_PAR,         ONLY : XUNDEF, NUNDEF
00064 USE MODD_ISBA_PAR,         ONLY : XWGMIN
00065 !
00066 USE MODI_GET_LUOUT
00067 USE MODI_UNPACK_SAME_RANK
00068 USE MODI_PACK_SAME_RANK
00069 USE MODI_ISBA_TO_TOPD
00070 USE MODI_RECHARGE_SURF_TOPD
00071 USE MODI_TOPODYN_LAT
00072 USE MODI_SAT_AREA_FRAC
00073 USE MODI_TOPD_TO_ISBA
00074 USE MODI_DIAG_ISBA_TO_ROUT
00075 USE MODI_ISBA_TO_TOPDSAT
00076 USE MODI_ROUTING
00077 USE MODI_OPEN_FILE
00078 USE MODI_WRITE_FILE_ISBAMAP
00079 USE MODI_CLOSE_FILE
00080 !
00081 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
00082 USE PARKIND1  ,ONLY : JPRB
00083 !
00084 IMPLICIT NONE
00085 !
00086 !*      0.1    declarations of arguments
00087 !
00088  CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes
00089  CHARACTER(LEN=*), INTENT(IN) :: HSTEP  ! atmospheric loop index
00090 INTEGER, INTENT(IN)          :: KI    ! Grid dimensions
00091 INTEGER, INTENT(IN)          :: KSTEP ! current time step 
00092 !
00093 !*      0.2    declarations of local variables
00094 !
00095 REAL, DIMENSION(SIZE(XRUNOFF_TOP,1)) :: ZWG2_TMP,ZWG3_TMP
00096 REAL, DIMENSION(SIZE(XRUNOFF_TOP,1)) :: ZWSAT_NAT
00097 REAL, DIMENSION(NNCAT,NMESHT) :: ZRT             ! recharge on TOP-LAT grid (m)
00098 REAL, DIMENSION(NNCAT,NMESHT) :: ZDEFT           ! local deficits on TOPODYN grid (m)
00099 REAL, DIMENSION(NNCAT,NMESHT) :: ZRI_WGIT        ! water changing of phase on TOPMODEL grid
00100 REAL, DIMENSION(NNCAT,NMESHT) :: ZRUNOFF_TOPD    ! Runoff on the Topodyn grid (m3/s)
00101 REAL, DIMENSION(NNCAT,NMESHT) :: ZDRAIN_TOPD     ! Drainage from Isba on Topodyn grid (m3/s)
00102 REAL, DIMENSION(NNCAT,NMESHT) :: ZKAPPA          ! topographic index
00103 REAL, DIMENSION(NNCAT)        :: ZKAPPAC         ! critical topographic index
00104 REAL, DIMENSION(KI)           :: ZRI             ! recharge on ISBA grid (m)
00105 REAL, DIMENSION(KI)           :: ZRI_WGI         ! water changing of phase on ISBA grid
00106 REAL, DIMENSION(KI)           :: Z_WSTOPI, Z_WFCTOPI
00107 REAL, DIMENSION(KI)           :: ZWM             ! Water content on SurfEx grid after the previous topodyn time step
00108 REAL, DIMENSION(KI)           :: ZWGIM           ! soil ice content at previous time step
00109 REAL, DIMENSION(KI)           :: ZRUNOFFC_FULL   ! Cumulated runoff from isba on the full domain (kg/m2)
00110 REAL, DIMENSION(KI)           :: ZRUNOFFC_FULLM  ! Cumulated runoff from isba on the full domain (kg/m2) at t-dt
00111 REAL, DIMENSION(KI)           :: ZRUNOFF_ISBA    ! Runoff from Isba (kg/m2)
00112 REAL, DIMENSION(KI)           :: ZDRAINC_FULL    ! Cumulated drainage from Isba on the full domain (kg/m2)
00113 REAL, DIMENSION(KI)           :: ZDRAINC_FULLM   ! Cumulated drainage from Isba on the full domain (kg/m2) at t-dt
00114 REAL, DIMENSION(KI)           :: ZDRAIN_ISBA     ! Drainage from Isba (m3/s)
00115 REAL, DIMENSION(KI)           :: ZDG_FULL
00116 REAL, DIMENSION(KI)           :: ZWG2_FULL, ZWG3_FULL, ZDG2_FULL, ZDG3_FULL
00117 REAL, DIMENSION(KI)           :: ZWGI_FULL
00118 REAL, DIMENSION(KI)           :: ZAS             ! Saturated area fraction for each Isba meshes
00119 REAL, DIMENSION(NNCAT)        :: Z_DW1,Z_DW2     ! Wsat-Wfc to actualise M in fonction of WI
00120 REAL                          :: ZAVG_MESH_SIZE
00121 LOGICAL, DIMENSION(NNCAT)     :: GTOPD           ! logical variable = true if topodyn_lat runs
00122 INTEGER                       :: JJ,JI           ! loop control 
00123 INTEGER                       :: IUNIT           ! unit number of results files
00124 INTEGER                       :: ILUOUT          ! unit number of listing file
00125 !
00126 REAL(KIND=JPRB) :: ZHOOK_HANDLE
00127 !-------------------------------------------------------------------------------
00128 IF (LHOOK) CALL DR_HOOK('COUPL_TOPD',0,ZHOOK_HANDLE)
00129 !
00130  CALL GET_LUOUT(HPROGRAM,ILUOUT)
00131 !
00132 !*       0.     Initialization:
00133 !               ---------------
00134 !
00135 ZWM(:) = XWG_FULL(:)
00136 ZWGIM(:) = XWGI_FULL(:)
00137 ZWG2_TMP(:) = XWG(:,2,1)
00138 ZWG3_TMP(:) = XWG(:,3,1)
00139 !
00140 !*       1.     ISBA => TOPODYN
00141 !               ---------------
00142 !*       1.1    Computation of the useful depth and water for lateral transfers
00143 !               -----------------------------------
00144 !
00145  CALL UNPACK_SAME_RANK(NR_NATURE,XDG(:,2,1),ZDG2_FULL)
00146  CALL UNPACK_SAME_RANK(NR_NATURE,XDG(:,3,1),ZDG3_FULL)
00147 WHERE ( ZDG2_FULL/=XUNDEF )
00148   ZDG_FULL = XFRAC_D2*ZDG2_FULL + XFRAC_D3*(ZDG3_FULL-ZDG2_FULL)
00149 ELSEWHERE
00150   ZDG_FULL = XUNDEF
00151 END WHERE
00152 !
00153  CALL UNPACK_SAME_RANK(NR_NATURE,XWG(:,2,1),ZWG2_FULL)
00154  CALL UNPACK_SAME_RANK(NR_NATURE,XWG(:,3,1),ZWG3_FULL)
00155 WHERE ( ZDG_FULL/=XUNDEF .AND. ZDG_FULL/=0. )
00156   XWG_FULL = XFRAC_D2*(ZDG2_FULL/ZDG_FULL)*ZWG2_FULL + XFRAC_D3*((ZDG3_FULL-ZDG2_FULL)/ZDG_FULL)*ZWG3_FULL
00157 ELSEWHERE
00158   XWG_FULL = XUNDEF
00159 END WHERE
00160 !
00161 !ludo prise en compte glace (pas de glace dans 3e couche)
00162  CALL UNPACK_SAME_RANK(NR_NATURE,XWGI(:,2,1),ZWGI_FULL)
00163 WHERE ( ZWGI_FULL/=XUNDEF .AND. XFRAC_D2>0 .AND. ZDG_FULL/=0. )
00164   XWGI_FULL = XFRAC_D2*(ZDG2_FULL/ZDG_FULL)*ZWGI_FULL
00165 ELSEWHERE
00166   XWGI_FULL = XUNDEF
00167 END WHERE
00168 !
00169 !
00170 WHERE ( XDTOPI/=XUNDEF ) 
00171   ZRI_WGI = ( XWGI_FULL - ZWGIM ) !m3/m3
00172 ELSEWHERE
00173   ZRI_WGI = 0.0
00174 END WHERE
00175 !
00176  CALL ISBA_TO_TOPD(ZRI_WGI,ZRI_WGIT)
00177 !
00178 !*       1.2    Water recharge 
00179 !               ---------------
00180 ! Topodyn uses :
00181 ! - a water recharge = water added since last time step to compute hydrological similarity indexes
00182 ! - the total water content to compute a deficit
00183 !
00184 ! This recharge is computed without regarding the changing of phase of water
00185 ! and the lateral transfers are performed regarding wsat et Wfc of last time step
00186 !
00187 WHERE ( XDTOPI/=XUNDEF )
00188   ZRI = ( (XWG_FULL - ZWM) + ZRI_WGI ) * XDTOPI
00189 ELSEWHERE
00190   ZRI = 0.0
00191 ENDWHERE
00192 !
00193 ! The water recharge on ISBA grid is computed on TOPMODEL grid
00194  CALL RECHARGE_SURF_TOPD(ZRI,ZRT,KI)
00195 !
00196 !*       2.     Lateral distribution
00197 !               --------------------
00198 !*       2.1    Computation of local deficits on TOPODYN grid
00199 !               ----------------------------------------
00200 !
00201  CALL TOPODYN_LAT(ZRT(:,:),ZDEFT(:,:),ZKAPPA(:,:),ZKAPPAC(:),GTOPD)
00202 !
00203 !*       2.2    Computation of contributive area on ISBA grid
00204 !               ----------------------------------------
00205 !
00206  CALL SAT_AREA_FRAC(ZDEFT,ZAS)
00207 !
00208  CALL PACK_SAME_RANK(NR_NATURE,ZAS,XAS_NATURE)
00209 !
00210 !*       3.    Deficit (m) -> water storage (m3/m3) and changing of phase
00211 !               ------------------------------------
00212 !ancien Wsat-Wfc pour actualisation param M
00213 Z_DW1(:)=0
00214 DO JJ=1,NNCAT
00215   Z_DW1(JJ) = SUM( XWSTOPT(JJ,:)-XWFCTOPT(JJ,:), MASK=XWSTOPT(JJ,:)/=XUNDEF ) / NNMC(JJ)
00216 ENDDO
00217 !
00218 DO JJ=1,NNCAT
00219   WHERE ( XDTOPT(JJ,:)/=XUNDEF .AND. XDTOPT(JJ,:)/=0. )
00220     XWTOPT(JJ,:) = XWSTOPT(JJ,:) - ( ZDEFT(JJ,:) / XDTOPT(JJ,:) )      
00221     !changing phase
00222     XWTOPT(JJ,:) = XWTOPT(JJ,:) - ZRI_WGIT(JJ,:)
00223   END WHERE
00224 ENDDO
00225 WHERE (XWTOPT > XWSTOPT ) XWTOPT = XWSTOPT
00226 !
00227 !actualisation Wsat, Wfc et Dmax pour prochain pas de temps
00228 WHERE ( ZWG2_FULL/=XUNDEF .AND. XWSTOPI/=0. )
00229   Z_WSTOPI  = XWSTOPI - XWGI_FULL
00230   Z_WFCTOPI = XWFCTOPI * Z_WSTOPI / XWSTOPI
00231 END WHERE
00232  CALL ISBA_TO_TOPD(Z_WSTOPI,XWSTOPT)
00233  CALL ISBA_TO_TOPD(Z_WFCTOPI,XWFCTOPT)
00234 !
00235 !ludo test empeche erreur num chgt phase
00236 WHERE ( ABS(XWSTOPT-XWTOPT) < 0.0000000001 ) XWSTOPT = XWTOPT
00237 !
00238 WHERE ( XWTOPT>XWSTOPT ) XWTOPT = XWSTOPT
00239 !
00240 WHERE ( XWFCTOPT/= XUNDEF ) XDMAXFC = (XWSTOPT - XWFCTOPT) * XDTOPT ! (m)
00241 !
00242 !actualisation M-> M=M*rapport Wsat-Wfc
00243 Z_DW2(:)=0
00244 DO JJ=1,NNCAT
00245   Z_DW2(JJ) = SUM( XWSTOPT(JJ,:)-XWFCTOPT(JJ,:), MASK=XWSTOPT(JJ,:)/=XUNDEF) / NNMC(JJ)
00246 ENDDO
00247 !
00248 XMPARA(:) = XMPARA(:) * Z_DW2(:)/Z_DW1(:)
00249 !
00250 !*       4.     TOPODYN => ISBA
00251 !               ---------------
00252 !*       4.1    Calculation of water storage on ISBA grid
00253 !               -----------------------------------------
00254 !
00255  CALL TOPD_TO_ISBA(KI,KSTEP,GTOPD)
00256  CALL PACK_SAME_RANK(NR_NATURE, (1-XFRAC_D2)*ZWG2_FULL + XFRAC_D2*XWG_FULL, XWG(:,2,1))
00257  CALL PACK_SAME_RANK(NR_NATURE, (1-XFRAC_D3)*ZWG3_FULL + XFRAC_D3*XWG_FULL, XWG(:,3,1))
00258 !
00259 !*       4.2    Budget correction
00260 !  -----------------------------------------
00261 !
00262 ZAVG_MESH_SIZE = SUM(XMESH_SIZE(:)) / SIZE(XMESH_SIZE,1)
00263 !
00264  CALL BUDG_WG(ZWG2_TMP(:),XWG(:,2,1),XDG(:,2,1),XMESH_SIZE,ZAVG_MESH_SIZE)
00265 !   
00266  CALL BUDG_WG(ZWG3_TMP(:),XWG(:,3,1),XDG(:,3,1)-XDG(:,2,1),XMESH_SIZE,ZAVG_MESH_SIZE)
00267 !
00268 !*      5.0    Total discharge
00269 !              ---------------
00270 !
00271 !*      5.1    Total water for runoff on TOPODYN grid
00272 !              ---------------------------------------
00273 !
00274  CALL PACK_SAME_RANK(NR_NATURE,XWSUPSAT,ZWSAT_NAT)
00275 WHERE ( ZWSAT_NAT(:)<10E-10 ) 
00276   ZWSAT_NAT(:) = 0.
00277 ELSEWHERE( ZWSAT_NAT(:)/=XUNDEF ) 
00278   XRUNOFF_TOP(:) = XRUNOFF_TOP(:) + ZWSAT_NAT(:)*XRHOLW*XDG(:,2,1)
00279 ENDWHERE
00280 !
00281 !write(*,*) ' ds coupl_topd ZWSAT_NAT',SUM(ZWSAT_NAT,MASK=ZWSAT_NAT(:)/=XUNDEF)!
00282  CALL UNPACK_SAME_RANK(NR_NATURE,XRUNOFF_TOP,ZRUNOFFC_FULL)
00283  CALL UNPACK_SAME_RANK(NR_NATURE,XAVG_RUNOFFCM,ZRUNOFFC_FULLM)
00284 !
00285  CALL DIAG_ISBA_TO_ROUT(ZRUNOFFC_FULL,ZRUNOFFC_FULLM,ZRUNOFF_ISBA)
00286 !
00287 XAVG_RUNOFFCM(:) = XRUNOFF_TOP(:)
00288 !
00289 WHERE (ZRUNOFF_ISBA==XUNDEF) ZRUNOFF_ISBA = 0.
00290 !
00291 ZRUNOFF_TOPD(:,:) = 0
00292 !
00293  CALL ISBA_TO_TOPDSAT(XKA_PRE,XKAC_PRE,KI,ZRUNOFF_ISBA,ZRUNOFF_TOPD)
00294 !
00295 !
00296 !*      5.2    Total water for drainage on TOPODYN grid
00297 !              ----------------------------------------
00298 !
00299  CALL UNPACK_SAME_RANK(NR_NATURE,XAVG_DRAINC*XATOP,ZDRAINC_FULL)
00300  CALL UNPACK_SAME_RANK(NR_NATURE,XAVG_DRAINCM*XATOP,ZDRAINC_FULLM)
00301 !
00302  CALL DIAG_ISBA_TO_ROUT(ZDRAINC_FULL,ZDRAINC_FULLM,ZDRAIN_ISBA)
00303 !
00304 WHERE (ZDRAIN_ISBA==XUNDEF) ZDRAIN_ISBA=0.
00305 !
00306 XAVG_DRAINCM(:)  = XAVG_DRAINC(:)
00307 !
00308 ZDRAIN_TOPD(:,:) = 0.0
00309 !
00310  CALL ISBA_TO_TOPD(ZDRAIN_ISBA,ZDRAIN_TOPD)
00311 !
00312 DO JJ=1,NNCAT
00313   DO JI=1,NNMC(JJ)
00314     IF (NMASKT(JJ,JI)/=NUNDEF) &
00315       ZDRAIN_TOPD(JJ,JI) = ZDRAIN_TOPD(JJ,JI) / NNPIX(NMASKT(JJ,JI))
00316   ENDDO
00317 ENDDO   
00318 !
00319 !*      6    Routing (runoff + drainage + exfiltration)
00320 !
00321  CALL ROUTING(ZRUNOFF_TOPD,ZDRAIN_TOPD,KSTEP)
00322 !
00323 XKA_PRE(:,:) = ZKAPPA(:,:)
00324 XKAC_PRE(:) = ZKAPPAC(:)
00325 !
00326 !*      7.0    Writing results in map files
00327 !              ----------------------------
00328 !
00329 IF (NFREQ_MAPS_ASAT/=0.AND.MOD(KSTEP,NFREQ_MAPS_ASAT)==0) THEN
00330   CALL OPEN_FILE('ASCII ',IUNIT,HFILE='carte_surfcont'//HSTEP,HFORM='FORMATTED',HACTION='WRITE')
00331   CALL WRITE_FILE_ISBAMAP(IUNIT,ZAS,KI)
00332   CALL CLOSE_FILE('ASCII ',IUNIT)
00333 ENDIF
00334 !
00335 IF (LHOOK) CALL DR_HOOK('COUPL_TOPD',1,ZHOOK_HANDLE)
00336 !
00337 CONTAINS
00338 !
00339 SUBROUTINE BUDG_WG(PWGM,PWG,PDG,PMESH_SIZE,PAVG_MESH_SIZE)
00340 !
00341 REAL, DIMENSION(:), INTENT(IN) :: PWGM
00342 REAL, DIMENSION(:), INTENT(INOUT) :: PWG
00343 REAL, DIMENSION(:), INTENT(IN) :: PDG
00344 REAL, DIMENSION(:), INTENT(IN) :: PMESH_SIZE
00345 REAL, INTENT(IN) :: PAVG_MESH_SIZE
00346 !
00347 REAL :: ZSTOCK_WGM, ZSTOCK_WG
00348 REAL :: ZAVG_DG, ZBUDG_WG
00349 REAL :: ZTMP, ZTMP2
00350 INTEGER :: JMESH
00351 !
00352 REAL(KIND=JPRB) :: ZHOOK_HANDLE
00353 !
00354 IF (LHOOK) CALL DR_HOOK('COUPL_TOPD:BUDG_WG',0,ZHOOK_HANDLE)
00355 !
00356 ZSTOCK_WGM = SUM( PWGM(:)*PDG(:)*PMESH_SIZE(:) )    ! water stocked in the ground (m3)
00357 ZSTOCK_WG  = SUM( PWG (:)*PDG(:)*PMESH_SIZE(:) )    ! water stocked in the layer2 (m3)
00358 ZAVG_DG = SUM(PDG(:)) / SIZE(PDG)
00359 ZBUDG_WG = ( ZSTOCK_WG - ZSTOCK_WGM )/ ZAVG_DG / PAVG_MESH_SIZE
00360 !ZTMP2=0.
00361 !DO JMESH=1,SIZE(PWG)             
00362 ! IF ((PWG(JMESH)/=XUNDEF).AND.(PWGM(JMESH)/=XUNDEF)) THEN
00363 !   IF(PWG(JMESH)/=PWGM(JMESH)) ZTMP2=ZTMP2+1.
00364 ! ENDIF
00365 !ENDDO
00366 !
00367 ZTMP  = COUNT( PWG(:)/=PWGM(:) )
00368 ZTMP2 = COUNT( PWG(:)/=PWGM(:) .AND. PWG(:)>XWGMIN+(ZBUDG_WG/ZTMP) )
00369 !   
00370 DO JMESH=1,SIZE(PWG)
00371   IF( PWG(JMESH)/=PWGM(JMESH) .AND. ZTMP2/=0. .AND. PWG(JMESH)>XWGMIN+(ZBUDG_WG/ZTMP) ) THEN
00372     !IF( PWG(JMESH)/=PWGM(JMESH) .AND. (ZTMP2/=0. ) THEN
00373     PWG(JMESH) = PWG(JMESH) - (ZBUDG_WG/ZTMP2)
00374     !PWG(JMESH) = MAX(PWG(JMESH) - (ZBUDG_WG/ZTMP2),XWGMIN)
00375   ENDIF
00376 ENDDO
00377 !write(*,*) KSTEP,COUNT(PWG(:)<XWGMIN),COUNT(PWG(:)<0.)
00378 !
00379 IF (LHOOK) CALL DR_HOOK('COUPL_TOPD:BUDG_WG',1,ZHOOK_HANDLE)
00380 !
00381 END SUBROUTINE BUDG_WG
00382 !
00383 END SUBROUTINE COUPL_TOPD