source: branches/publications/ORCHIDEE-PEAT_r5488/src_sechiba/hydro_subgrid.f90 @ 8787

Last change on this file since 8787 was 3342, checked in by albert.jornet, 9 years ago

New+fix: TOPMODEL debug and TOPMODEL pass exam
File size: 9.4 KB
Line 
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!     #################
30MODULE 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
40CONTAINS
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!
55INTEGER(i_std), INTENT(in)                         :: kjpindex         !! Domain size
56REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ruu_ch
57REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: PM
58REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: PTI_MIN
59REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: PTI_MAX
60REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: PAS
61REAL(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
67REAL(r_std), INTENT (in)       :: PD_TOP
68
69!
70REAL(r_std), DIMENSION(kjpindex), INTENT(OUT)  :: PFSAT,PFWET,PWT1,PWT2,PWT3,PWT4
71!                                   PFSAT   = satured fraction
72!                                   PFWET   = wetland fraction
73!
74REAL(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!
79REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: phumtot
80REAL(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!
85REAL, 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!
89REAL, DIMENSION(:,:), ALLOCATABLE     :: ZCOMP,ZCOMP_WET
90!                                                      ZCOMP     = array for wt comparaison with wt_array
91!
92INTEGER, 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!
95INTEGER                               :: INI, IND, I 
96!                                        INI      = nombre de mailles       
97!                                        IND      = Indice de boucle       
98!REAL(r_std), DIMENSION(kjpindex)                     :: PW_TOP, PWI_TOP
99REAL(r_std), DIMENSION (kjpindex)                    :: phumtot_ave
100INTEGER(i_std)                                      :: ji,jv
101REAL(r_std), DIMENSION (kjpindex)                    :: pfrac_froz_fsat
102INTEGER(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!
118INI=kjpindex
119ALLOCATE(ZCOMP(INI,SIZE(PTAB_WTOP,2)))
120ALLOCATE(ZCOMP_WET(INI,SIZE(PTAB_WTOP_WET,2)))
121!
122I_IND (:) =0 !indice seuil pr calcul de fraction a saturation
123I_IND_WET (:) =0 !indice seuil pr fraction de wetland (avc WTD=0)
124I_IND_WT1 (:) =0 !indice des fractions de wetland avec une WTD en-dessous du sol
125I_IND_WT2 (:) =0
126I_IND_WT3 (:) =0
127I_IND_WT4 (:) =0
128
129ZW_TOP(:) =0.0
130ZW_TOP_WET(:) =0.0
131PFSAT (:) =0.0
132PFWET (:) =0.0
133!
134!       
135IND=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))
141ZW_TOP(:) = MIN((phumtot_ave(:))/(PD_TOP*1000.), PTAB_WTOP(:,1))
142ZW_TOP_WET(:) = MIN((phumtot_ave(:))/(PD_TOP*1000.), PTAB_WTOP_WET(:,1))
143
144!
145! compare wt_array and WT
146!
147ZCOMP(:,:) = 999.0
148ZCOMP_WET(:,:) = 999.0
149!
150DO 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(:))
158ENDDO
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
163I_IND(:)=INT(MINLOC(ZCOMP(:,:),2))
164I_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
172WHERE ((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
181elsewhere
182    I_IND_WT1(:)=999
183    I_IND_WT2(:)=999
184    I_IND_WT3(:)=999
185    I_IND_WT4(:)=999
186endwhere
187
188WHERE ( I_IND_WT1(:) .LE. 0.0 ) 
189    I_IND_WT1(:)=1
190endwhere
191WHERE ( I_IND_WT2(:) .LE. 0.0 ) 
192    I_IND_WT2(:)=1
193endwhere
194WHERE ( I_IND_WT3(:) .LE. 0.0 ) 
195    I_IND_WT3(:)=1
196endwhere
197WHERE ( I_IND_WT4(:) .LE. 0.0 ) 
198    I_IND_WT4(:)=1
199endwhere
200
201DO I=1,INI
202  nb_count(I)=COUNT(PTAB_FSAT (I,:)>0.0)!!= COUNT(PTAB_FWET (I,:)>0.0) aussi
203ENDDO
204
205DO 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
223ENDDO
224!
225
226WHERE ( 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
233endwhere
234
235!prise en compte du gel pour diminuer les differentes fractions
236WHERE ( 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(:)
243ENDWHERE
244
245!
246!-------------------------------------------------------------------------------
247!
248END SUBROUTINE hydro_subgrid_main
249END MODULE hydro_subgrid 
Note: See TracBrowser for help on using the repository browser.