| 1 | ! ================================================================================================================================= |
|---|
| 2 | ! MODULE : hydro_subgrid |
|---|
| 3 | ! |
|---|
| 4 | ! CONTACT : orchidee-help _at_ ipsl.jussieu.fr |
|---|
| 5 | ! |
|---|
| 6 | ! LICENCE : IPSL (2006) |
|---|
| 7 | ! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC |
|---|
| 8 | ! |
|---|
| 9 | !>\BRIEF This module computes TOPMODEL saturated fraction. |
|---|
| 10 | !! |
|---|
| 11 | !!\n DESCRIPTION : contains. |
|---|
| 12 | !! |
|---|
| 13 | !! RECENT CHANGE(S) : None |
|---|
| 14 | !! |
|---|
| 15 | !! REFERENCE(S) : None |
|---|
| 16 | !! |
|---|
| 17 | !! SVN : |
|---|
| 18 | !! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/perso/shushi.peng/ORCHIDEE/src_sechiba/hydro_subgrid.f90 $ |
|---|
| 19 | !! $Date: 2015-03-10 10:16:04 +0100 (Tue, 10 Mar 2015) $ |
|---|
| 20 | !! $Revision: 2535 $ |
|---|
| 21 | !! \n |
|---|
| 22 | !_ ================================================================================================================================ |
|---|
| 23 | |
|---|
| 24 | !----------------------------------------------------------------- |
|---|
| 25 | !--------------- special set of characters for SCCS information |
|---|
| 26 | !----------------------------------------------------------------- |
|---|
| 27 | ! %Z% Lib:%F%, Version:%I%, Date:%D%, Last modified:%E% |
|---|
| 28 | !----------------------------------------------------------------- |
|---|
| 29 | ! ################# |
|---|
| 30 | MODULE hydro_subgrid |
|---|
| 31 | ! ################# |
|---|
| 32 | ! |
|---|
| 33 | USE ioipsl |
|---|
| 34 | USE constantes |
|---|
| 35 | USE constantes_soil |
|---|
| 36 | !pss! USE constantes_veg |
|---|
| 37 | !! USE reqdprec |
|---|
| 38 | ! |
|---|
| 39 | IMPLICIT NONE |
|---|
| 40 | CONTAINS |
|---|
| 41 | ! |
|---|
| 42 | SUBROUTINE hydro_subgrid_main(kjpindex, PTAB_FSAT, PTAB_WTOP, phumtot, profil_froz_hydro, PFSAT,& |
|---|
| 43 | & PTAB_FWET,PTAB_WTOP_WET,PFWET, PD_TOP, & |
|---|
| 44 | & ruu_ch, PWT1,PWT2,PWT3,PWT4,PM,PTI_MIN,PTI_MAX,PAS, dz) |
|---|
| 45 | ! |
|---|
| 46 | ! |
|---|
| 47 | !------------------------------------------------------------------------------- |
|---|
| 48 | ! |
|---|
| 49 | !* 0. DECLARATIONS |
|---|
| 50 | ! ------------ |
|---|
| 51 | ! |
|---|
| 52 | !* 0.1 declarations of arguments |
|---|
| 53 | ! |
|---|
| 54 | ! |
|---|
| 55 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size |
|---|
| 56 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ruu_ch |
|---|
| 57 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: PM |
|---|
| 58 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: PTI_MIN |
|---|
| 59 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: PTI_MAX |
|---|
| 60 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: PAS |
|---|
| 61 | REAL(r_std),DIMENSION (nslm+1), INTENT (in) :: dz |
|---|
| 62 | ! PW_TOP = total liquid volumetric water |
|---|
| 63 | ! content of the upper layer |
|---|
| 64 | ! PWI_TOP = total ice volumetric water |
|---|
| 65 | ! content of the below layer |
|---|
| 66 | |
|---|
| 67 | REAL(r_std), INTENT (in) :: PD_TOP |
|---|
| 68 | |
|---|
| 69 | ! |
|---|
| 70 | REAL(r_std), DIMENSION(kjpindex), INTENT(OUT) :: PFSAT,PFWET,PWT1,PWT2,PWT3,PWT4 |
|---|
| 71 | ! PFSAT = satured fraction |
|---|
| 72 | ! PFWET = wetland fraction |
|---|
| 73 | ! |
|---|
| 74 | REAL(r_std), DIMENSION(kjpindex,1000), INTENT(IN) :: PTAB_FSAT, PTAB_WTOP, PTAB_FWET, PTAB_WTOP_WET |
|---|
| 75 | ! PTAB_FSAT = Satured fraction array |
|---|
| 76 | ! PTAB_WTOP = Active TOPMODEL-layer array |
|---|
| 77 | ! PTAB_FWET = Wetland fraction array |
|---|
| 78 | ! |
|---|
| 79 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: phumtot |
|---|
| 80 | REAL(r_std),DIMENSION (kjpindex, nslm), INTENT (in) :: profil_froz_hydro |
|---|
| 81 | !REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in) :: psoiltile |
|---|
| 82 | |
|---|
| 83 | !* 0.2 declarations of local variables |
|---|
| 84 | ! |
|---|
| 85 | REAL, DIMENSION(kjpindex) :: ZW_TOP,ZW_TOP_WET |
|---|
| 86 | !REAL, DIMENSION(kjpindex) :: ZFSAT, ZW_TOP, ZFWET |
|---|
| 87 | ! ZW_TOP = ative TOPMODEL-soil moisture at 't' (m3 m-3) |
|---|
| 88 | ! |
|---|
| 89 | REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOMP,ZCOMP_WET |
|---|
| 90 | ! ZCOMP = array for wt comparaison with wt_array |
|---|
| 91 | ! |
|---|
| 92 | INTEGER, DIMENSION(kjpindex) :: I_IND,I_IND_WET,I_IND_WT1,I_IND_WT2,I_IND_WT3,I_IND_WT4 |
|---|
| 93 | ! I_IND = change in xsat (or fsat) index |
|---|
| 94 | ! |
|---|
| 95 | INTEGER :: INI, IND, I |
|---|
| 96 | ! INI = nombre de mailles |
|---|
| 97 | ! IND = Indice de boucle |
|---|
| 98 | !REAL(r_std), DIMENSION(kjpindex) :: PW_TOP, PWI_TOP |
|---|
| 99 | REAL(r_std), DIMENSION (kjpindex) :: phumtot_ave |
|---|
| 100 | INTEGER(i_std) :: ji,jv |
|---|
| 101 | REAL(r_std), DIMENSION (kjpindex) :: pfrac_froz_fsat |
|---|
| 102 | INTEGER(r_std), DIMENSION (kjpindex) :: nb_count |
|---|
| 103 | !------------------------------------------------------------------------------- |
|---|
| 104 | ! |
|---|
| 105 | !valeur pr chaque grid-cell |
|---|
| 106 | phumtot_ave(:) = 0.0 |
|---|
| 107 | ! |
|---|
| 108 | DO ji = 1, kjpindex |
|---|
| 109 | !!pmc_ave(ji) = SUM(pmc(ji,:,1)*psoiltile(ji,1)+pmc(ji,:,2)*psoiltile(ji,2)+pmc(ji,:,3)*psoiltile(ji,3))!SUM sur les 11 couches |
|---|
| 110 | !variable _ave servait initialement a faire la somme sur les PFT ou sur les types de sol. inutile maintenant |
|---|
| 111 | phumtot_ave(ji) = phumtot(ji) |
|---|
| 112 | ENDDO |
|---|
| 113 | |
|---|
| 114 | DO ji = 1, kjpindex |
|---|
| 115 | pfrac_froz_fsat(ji)= (SUM(profil_froz_hydro(ji,1:9) * dz(1:9)) / SUM(dz(1:9))) |
|---|
| 116 | ENDDO |
|---|
| 117 | ! |
|---|
| 118 | INI=kjpindex |
|---|
| 119 | ALLOCATE(ZCOMP(INI,SIZE(PTAB_WTOP,2))) |
|---|
| 120 | ALLOCATE(ZCOMP_WET(INI,SIZE(PTAB_WTOP_WET,2))) |
|---|
| 121 | ! |
|---|
| 122 | I_IND (:) =0 !indice seuil pr calcul de fraction a saturation |
|---|
| 123 | I_IND_WET (:) =0 !indice seuil pr fraction de wetland (avc WTD=0) |
|---|
| 124 | I_IND_WT1 (:) =0 !indice des fractions de wetland avec une WTD en-dessous du sol |
|---|
| 125 | I_IND_WT2 (:) =0 |
|---|
| 126 | I_IND_WT3 (:) =0 |
|---|
| 127 | I_IND_WT4 (:) =0 |
|---|
| 128 | |
|---|
| 129 | ZW_TOP(:) =0.0 |
|---|
| 130 | ZW_TOP_WET(:) =0.0 |
|---|
| 131 | PFSAT (:) =0.0 |
|---|
| 132 | PFWET (:) =0.0 |
|---|
| 133 | ! |
|---|
| 134 | ! |
|---|
| 135 | IND=SIZE(PTAB_WTOP,2) |
|---|
| 136 | ! |
|---|
| 137 | ! Here, we determinate the new fsat by comparaison between wt_array and WT |
|---|
| 138 | ! |
|---|
| 139 | !ZW_TOP(:) = MIN(((PW_TOP(:) + PWI_TOP(:)+pdrainage_t(:)))/(PD_TOP(:)*1000.), PTAB_WTOP(:,1)) |
|---|
| 140 | !ZW_TOP_WET(:) = MIN(((PW_TOP(:) + PWI_TOP(:)+pdrainage_t(:)))/(PD_TOP(:)*1000.), PTAB_WTOP_WET(:,1)) |
|---|
| 141 | ZW_TOP(:) = MIN((phumtot_ave(:))/(PD_TOP*1000.), PTAB_WTOP(:,1)) |
|---|
| 142 | ZW_TOP_WET(:) = MIN((phumtot_ave(:))/(PD_TOP*1000.), PTAB_WTOP_WET(:,1)) |
|---|
| 143 | |
|---|
| 144 | ! |
|---|
| 145 | ! compare wt_array and WT |
|---|
| 146 | ! |
|---|
| 147 | ZCOMP(:,:) = 999.0 |
|---|
| 148 | ZCOMP_WET(:,:) = 999.0 |
|---|
| 149 | ! |
|---|
| 150 | DO I=1,IND |
|---|
| 151 | !prise en compte du gel pour le calcul de fsat: |
|---|
| 152 | !soluce 1: diminue la quantité d eau utilise pour calculer le deficit |
|---|
| 153 | !soluce 2: utilise la qte d eau totale (liq + solide) puis diminue le fsat obtenu avec ce calcul |
|---|
| 154 | ! ZCOMP(:,I) = ABS(PTAB_WTOP(:,I) - ((pfrac_froz_gqsb(:)* dsg_ave(:) + & |
|---|
| 155 | ! & (PD_TOP(:)-dsg_ave(:))*pfrac_froz_bqsb(:))/PD_TOP(:))*(ruu_ch(:)/1000.0) - ZW_TOP(:)) |
|---|
| 156 | ZCOMP(:,I) = ABS(PTAB_WTOP(:,I) - ZW_TOP(:)) |
|---|
| 157 | ZCOMP_WET(:,I) = ABS(PTAB_WTOP_WET(:,I) - ZW_TOP_WET(:)) |
|---|
| 158 | ENDDO |
|---|
| 159 | ! |
|---|
| 160 | ! calculate array index where tab_wt = WT |
|---|
| 161 | ! |
|---|
| 162 | !recherche de l indice seuil pour les fractions a saturation et fractions de wetland |
|---|
| 163 | I_IND(:)=INT(MINLOC(ZCOMP(:,:),2)) |
|---|
| 164 | I_IND_WET(:)=INT(MINLOC(ZCOMP_WET(:,:),2)) |
|---|
| 165 | |
|---|
| 166 | !boucle de type ci-dessus peut etre utilise si la fonction MINLOC n est pas vectorise par le compilo |
|---|
| 167 | !DO ji = 1, kjpindex |
|---|
| 168 | ! I_IND(ji)=INT(MINLOC(ZCOMP(ji,:),1)) |
|---|
| 169 | !ENDDO |
|---|
| 170 | |
|---|
| 171 | |
|---|
| 172 | WHERE ((PTI_MAX(:) /= -99.99 ) .AND. & |
|---|
| 173 | & (ruu_ch(:) .GE. min_sechiba) .AND. (PD_TOP .GE. min_sechiba) .AND. (PAS(:) .GE. min_sechiba)) |
|---|
| 174 | I_IND_WT1(:)=INT(I_IND_WET(:) - ((4*WTD1_borne/PD_TOP)/(ruu_ch(:)*PD_TOP/1000.))/PAS(:))! indice pour deficit entre 0 et -6cm |
|---|
| 175 | I_IND_WT2(:)=INT(I_IND_WET(:) - ((4*WTD2_borne/PD_TOP)/(ruu_ch(:)*PD_TOP/1000.))/PAS(:))! indice pour deficit entre 0 et -12 cm |
|---|
| 176 | I_IND_WT3(:)=INT(I_IND_WET(:) - ((4*WTD3_borne/PD_TOP)/(ruu_ch(:)*PD_TOP/1000.))/PAS(:))! indice pour deficit entre 0 et -18cm |
|---|
| 177 | I_IND_WT4(:)=INT(I_IND_WET(:) - ((4*WTD4_borne/PD_TOP)/(ruu_ch(:)*PD_TOP/1000.))/PAS(:))! indice pour deficit entre 0 et -24cm |
|---|
| 178 | ! le pd_top au denominateur est un peu etrange |
|---|
| 179 | ! pourtant ce calcul donne les bonnes fractions: |
|---|
| 180 | ! un test a ete realise en utilisant I_IND et en calculant la fraction pour un deficit entre 0 et 2 m => fraction ~ 1 |
|---|
| 181 | elsewhere |
|---|
| 182 | I_IND_WT1(:)=999 |
|---|
| 183 | I_IND_WT2(:)=999 |
|---|
| 184 | I_IND_WT3(:)=999 |
|---|
| 185 | I_IND_WT4(:)=999 |
|---|
| 186 | endwhere |
|---|
| 187 | |
|---|
| 188 | WHERE ( I_IND_WT1(:) .LE. 0.0 ) |
|---|
| 189 | I_IND_WT1(:)=1 |
|---|
| 190 | endwhere |
|---|
| 191 | WHERE ( I_IND_WT2(:) .LE. 0.0 ) |
|---|
| 192 | I_IND_WT2(:)=1 |
|---|
| 193 | endwhere |
|---|
| 194 | WHERE ( I_IND_WT3(:) .LE. 0.0 ) |
|---|
| 195 | I_IND_WT3(:)=1 |
|---|
| 196 | endwhere |
|---|
| 197 | WHERE ( I_IND_WT4(:) .LE. 0.0 ) |
|---|
| 198 | I_IND_WT4(:)=1 |
|---|
| 199 | endwhere |
|---|
| 200 | |
|---|
| 201 | DO I=1,INI |
|---|
| 202 | nb_count(I)=COUNT(PTAB_FSAT (I,:)>0.0)!!= COUNT(PTAB_FWET (I,:)>0.0) aussi |
|---|
| 203 | ENDDO |
|---|
| 204 | |
|---|
| 205 | DO I=1,INI |
|---|
| 206 | ! calculate fsat |
|---|
| 207 | IF(nb_count(I)>0.0)THEN |
|---|
| 208 | PFSAT(I) = PTAB_FSAT (I,I_IND(I)) |
|---|
| 209 | PFWET(I) = PTAB_FWET (I,I_IND_WET(I)) |
|---|
| 210 | !PFWET(I) = PTAB_FWET (I,I_IND(I)) ! test pour voir si somme des fractions de deficit entre 0 et 2m ~ 1 |
|---|
| 211 | PWT1(I) = PTAB_FWET (I,I_IND_WT1(I))-PFWET(I) |
|---|
| 212 | PWT2(I) = PTAB_FWET (I,I_IND_WT2(I))-(PWT1(I)+PFWET(I)) |
|---|
| 213 | PWT3(I) = PTAB_FWET (I,I_IND_WT3(I))-(PWT2(I)+PWT1(I)+PFWET(I)) |
|---|
| 214 | PWT4(I) = PTAB_FWET (I,I_IND_WT4(I))-(PWT3(I)+PWT2(I)+PWT1(I)+PFWET(I)) |
|---|
| 215 | ELSE |
|---|
| 216 | PFSAT(I) = 0.0 !Australia |
|---|
| 217 | PFWET(I) = 0.0 |
|---|
| 218 | PWT1(I) = 0.0 |
|---|
| 219 | PWT2(I) = 0.0 |
|---|
| 220 | PWT3(I) = 0.0 |
|---|
| 221 | PWT4(I) = 0.0 |
|---|
| 222 | ENDIF |
|---|
| 223 | ENDDO |
|---|
| 224 | ! |
|---|
| 225 | |
|---|
| 226 | WHERE ( PTAB_WTOP(:,1) .EQ. PTAB_WTOP(:,900) ) |
|---|
| 227 | PFSAT(:)=0.0 |
|---|
| 228 | PFWET(:)=0.0 |
|---|
| 229 | PWT1(:)=0.0 |
|---|
| 230 | PWT2(:)=0.0 |
|---|
| 231 | PWT3(:)=0.0 |
|---|
| 232 | PWT4(:)=0.0 |
|---|
| 233 | endwhere |
|---|
| 234 | |
|---|
| 235 | !prise en compte du gel pour diminuer les differentes fractions |
|---|
| 236 | WHERE ( phumtot_ave(:).GT. 0.0 ) |
|---|
| 237 | PFSAT(:)=(1-pfrac_froz_fsat(:))*PFSAT(:) |
|---|
| 238 | PFWET(:)=(1-pfrac_froz_fsat(:))*PFWET(:) |
|---|
| 239 | PWT1(:)= (1-pfrac_froz_fsat(:))*PWT1(:) |
|---|
| 240 | PWT2(:)= (1-pfrac_froz_fsat(:))*PWT2(:) |
|---|
| 241 | PWT3(:)= (1-pfrac_froz_fsat(:))*PWT3(:) |
|---|
| 242 | PWT4(:)= (1-pfrac_froz_fsat(:))*PWT4(:) |
|---|
| 243 | ENDWHERE |
|---|
| 244 | |
|---|
| 245 | ! |
|---|
| 246 | !------------------------------------------------------------------------------- |
|---|
| 247 | ! |
|---|
| 248 | END SUBROUTINE hydro_subgrid_main |
|---|
| 249 | END MODULE hydro_subgrid |
|---|
