source: branches/publications/ORCHIDEE-MICT-OP-r6850/src_sechiba/topmodel_subgrid.f90 @ 7421

Last change on this file since 7421 was 5617, checked in by albert.jornet, 6 years ago

Refactor: moved missing TopModel? declarations (IOISPL history and its parameters )
File size: 9.4 KB
Line 
1! =================================================================================================================================
2! MODULE        : topmodel_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 topmodel_subgrid 
31!     #################
32!
33      USE ioipsl_para
34      USE constantes_soil
35      USE topmodel_var
36!
37      IMPLICIT NONE
38CONTAINS
39!
40      SUBROUTINE topmodel_subgrid_main(kjpindex, PTAB_FSAT, PTAB_WTOP, phumtot, profil_froz_hydro, PFSAT,&
41    & PTAB_FWET,PTAB_WTOP_WET,PFWET, PD_TOP, &
42    & ruu_ch, PWT1,PWT2,PWT3,PWT4,PM,PTI_MIN,PTI_MAX,PAS, dz)
43!
44!
45!-------------------------------------------------------------------------------
46!
47!*       0.     DECLARATIONS
48!               ------------
49!
50!*      0.1    declarations of arguments
51!
52!
53INTEGER(i_std), INTENT(in)                         :: kjpindex         !! Domain size
54REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ruu_ch
55REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: PM
56REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: PTI_MIN
57REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: PTI_MAX
58REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: PAS
59REAL(r_std),DIMENSION (nslm+1),   INTENT (in)      :: dz
60!                                   PW_TOP  = total liquid volumetric water
61!                                             content of the upper layer
62!                                   PWI_TOP = total ice volumetric water
63!                                             content of the below layer
64
65REAL(r_std), INTENT (in)       :: PD_TOP
66
67!
68REAL(r_std), DIMENSION(kjpindex), INTENT(OUT)  :: PFSAT,PFWET,PWT1,PWT2,PWT3,PWT4
69!                                   PFSAT   = satured fraction
70!                                   PFWET   = wetland fraction
71!
72REAL(r_std), DIMENSION(kjpindex,1000), INTENT(IN) :: PTAB_FSAT, PTAB_WTOP, PTAB_FWET, PTAB_WTOP_WET 
73!                                   PTAB_FSAT = Satured fraction array
74!                                   PTAB_WTOP = Active TOPMODEL-layer array
75!                                   PTAB_FWET = Wetland fraction array
76!
77REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: phumtot
78REAL(r_std),DIMENSION (kjpindex, nslm), INTENT (in) :: profil_froz_hydro 
79!REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in) :: psoiltile
80
81!*      0.2    declarations of local variables
82!
83REAL, DIMENSION(kjpindex)                     :: ZW_TOP,ZW_TOP_WET
84!REAL, DIMENSION(kjpindex)                     :: ZFSAT, ZW_TOP, ZFWET
85!                                        ZW_TOP  = ative TOPMODEL-soil moisture at 't' (m3 m-3)
86!
87REAL, DIMENSION(:,:), ALLOCATABLE     :: ZCOMP,ZCOMP_WET
88!                                                      ZCOMP     = array for wt comparaison with wt_array
89!
90INTEGER, DIMENSION(kjpindex)                  :: I_IND,I_IND_WET,I_IND_WT1,I_IND_WT2,I_IND_WT3,I_IND_WT4
91!                                        I_IND = change in xsat (or fsat) index
92!
93INTEGER                               :: INI, IND, I 
94!                                        INI      = nombre de mailles       
95!                                        IND      = Indice de boucle       
96!REAL(r_std), DIMENSION(kjpindex)                     :: PW_TOP, PWI_TOP
97REAL(r_std), DIMENSION (kjpindex)                    :: phumtot_ave
98INTEGER(i_std)                                      :: ji,jv
99REAL(r_std), DIMENSION (kjpindex)                    :: pfrac_froz_fsat
100INTEGER(r_std), DIMENSION (kjpindex)                 :: nb_count
101!-------------------------------------------------------------------------------
102!
103!valeur pr chaque grid-cell
104    phumtot_ave(:) = 0.0
105!
106    DO ji = 1, kjpindex
107       !!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
108       !variable _ave servait initialement a faire la somme sur les PFT ou sur les types de sol. inutile maintenant
109       phumtot_ave(ji) = phumtot(ji)     
110    ENDDO
111
112    DO ji = 1, kjpindex
113       pfrac_froz_fsat(ji)= (SUM(profil_froz_hydro(ji,1:9) * dz(1:9)) / SUM(dz(1:9)))
114    ENDDO
115!
116INI=kjpindex
117ALLOCATE(ZCOMP(INI,SIZE(PTAB_WTOP,2)))
118ALLOCATE(ZCOMP_WET(INI,SIZE(PTAB_WTOP_WET,2)))
119!
120I_IND (:) =0 !indice seuil pr calcul de fraction a saturation
121I_IND_WET (:) =0 !indice seuil pr fraction de wetland (avc WTD=0)
122I_IND_WT1 (:) =0 !indice des fractions de wetland avec une WTD en-dessous du sol
123I_IND_WT2 (:) =0
124I_IND_WT3 (:) =0
125I_IND_WT4 (:) =0
126
127ZW_TOP(:) =0.0
128ZW_TOP_WET(:) =0.0
129PFSAT (:) =0.0
130PFWET (:) =0.0
131!
132!       
133IND=SIZE(PTAB_WTOP,2)
134! 
135! Here, we determinate the new fsat by comparaison between wt_array and WT
136!
137!ZW_TOP(:) = MIN(((PW_TOP(:) + PWI_TOP(:)+pdrainage_t(:)))/(PD_TOP(:)*1000.), PTAB_WTOP(:,1))
138!ZW_TOP_WET(:) = MIN(((PW_TOP(:) + PWI_TOP(:)+pdrainage_t(:)))/(PD_TOP(:)*1000.), PTAB_WTOP_WET(:,1))
139ZW_TOP(:) = MIN((phumtot_ave(:))/(PD_TOP*1000.), PTAB_WTOP(:,1))
140ZW_TOP_WET(:) = MIN((phumtot_ave(:))/(PD_TOP*1000.), PTAB_WTOP_WET(:,1))
141
142!
143! compare wt_array and WT
144!
145ZCOMP(:,:) = 999.0
146ZCOMP_WET(:,:) = 999.0
147!
148DO I=1,IND
149  !prise en compte du gel pour le calcul de fsat:
150  !soluce 1: diminue la quantité d eau utilise pour calculer le deficit
151  !soluce 2: utilise la qte d eau totale (liq + solide) puis diminue le fsat obtenu avec ce calcul
152  ! ZCOMP(:,I) = ABS(PTAB_WTOP(:,I) - ((pfrac_froz_gqsb(:)* dsg_ave(:) + &
153  !    & (PD_TOP(:)-dsg_ave(:))*pfrac_froz_bqsb(:))/PD_TOP(:))*(ruu_ch(:)/1000.0) - ZW_TOP(:))
154  ZCOMP(:,I) = ABS(PTAB_WTOP(:,I) - ZW_TOP(:))
155  ZCOMP_WET(:,I) = ABS(PTAB_WTOP_WET(:,I) - ZW_TOP_WET(:))
156ENDDO
157!
158! calculate array index where tab_wt = WT
159!
160!recherche de l indice seuil pour les fractions a saturation et fractions de wetland
161I_IND(:)=INT(MINLOC(ZCOMP(:,:),2))
162I_IND_WET(:)=INT(MINLOC(ZCOMP_WET(:,:),2))
163
164!boucle de type ci-dessus peut etre utilise si la fonction MINLOC n est pas vectorise par le compilo
165!DO ji = 1, kjpindex
166!   I_IND(ji)=INT(MINLOC(ZCOMP(ji,:),1))
167!ENDDO
168
169
170WHERE ((PTI_MAX(:) /= -99.99 ) .AND. &
171        & (ruu_ch(:) .GE. min_sechiba) .AND. (PD_TOP .GE. min_sechiba) .AND. (PAS(:) .GE. min_sechiba))
172    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
173    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
174    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
175    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
176    ! le pd_top au denominateur est un peu etrange
177    ! pourtant ce calcul donne les bonnes fractions:
178    ! un test a ete realise en utilisant I_IND et en calculant la fraction pour un deficit entre 0 et 2 m => fraction ~ 1
179elsewhere
180    I_IND_WT1(:)=999
181    I_IND_WT2(:)=999
182    I_IND_WT3(:)=999
183    I_IND_WT4(:)=999
184endwhere
185
186WHERE ( I_IND_WT1(:) .LE. 0.0 ) 
187    I_IND_WT1(:)=1
188endwhere
189WHERE ( I_IND_WT2(:) .LE. 0.0 ) 
190    I_IND_WT2(:)=1
191endwhere
192WHERE ( I_IND_WT3(:) .LE. 0.0 ) 
193    I_IND_WT3(:)=1
194endwhere
195WHERE ( I_IND_WT4(:) .LE. 0.0 ) 
196    I_IND_WT4(:)=1
197endwhere
198
199DO I=1,INI
200  nb_count(I)=COUNT(PTAB_FSAT (I,:)>0.0)!!= COUNT(PTAB_FWET (I,:)>0.0) aussi
201ENDDO
202
203DO I=1,INI
204!  calculate fsat
205  IF(nb_count(I)>0.0)THEN
206      PFSAT(I) = PTAB_FSAT (I,I_IND(I))
207      PFWET(I) = PTAB_FWET (I,I_IND_WET(I))
208      !PFWET(I) = PTAB_FWET (I,I_IND(I)) ! test pour voir si somme des fractions de deficit entre 0 et 2m ~ 1
209      PWT1(I) = PTAB_FWET (I,I_IND_WT1(I))-PFWET(I) 
210      PWT2(I) = PTAB_FWET (I,I_IND_WT2(I))-(PWT1(I)+PFWET(I))
211      PWT3(I) = PTAB_FWET (I,I_IND_WT3(I))-(PWT2(I)+PWT1(I)+PFWET(I))
212      PWT4(I) = PTAB_FWET (I,I_IND_WT4(I))-(PWT3(I)+PWT2(I)+PWT1(I)+PFWET(I))
213   ELSE
214       PFSAT(I) = 0.0 !Australia
215       PFWET(I) = 0.0
216       PWT1(I) = 0.0
217       PWT2(I) = 0.0
218       PWT3(I) = 0.0
219       PWT4(I) = 0.0
220   ENDIF
221ENDDO
222!
223
224WHERE ( PTAB_WTOP(:,1) .EQ. PTAB_WTOP(:,900) ) 
225    PFSAT(:)=0.0
226    PFWET(:)=0.0
227    PWT1(:)=0.0
228    PWT2(:)=0.0
229    PWT3(:)=0.0
230    PWT4(:)=0.0
231endwhere
232
233!prise en compte du gel pour diminuer les differentes fractions
234WHERE ( phumtot_ave(:).GT. 0.0 )
235    PFSAT(:)=(1-pfrac_froz_fsat(:))*PFSAT(:)
236    PFWET(:)=(1-pfrac_froz_fsat(:))*PFWET(:) 
237    PWT1(:)= (1-pfrac_froz_fsat(:))*PWT1(:)
238    PWT2(:)= (1-pfrac_froz_fsat(:))*PWT2(:)
239    PWT3(:)= (1-pfrac_froz_fsat(:))*PWT3(:)
240    PWT4(:)= (1-pfrac_froz_fsat(:))*PWT4(:)
241ENDWHERE
242
243!
244!-------------------------------------------------------------------------------
245!
246END SUBROUTINE topmodel_subgrid_main
247END MODULE topmodel_subgrid 
Note: See TracBrowser for help on using the repository browser.