/[lmdze]/trunk/Sources/phylmd/clmain.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/clmain.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/phylmd/clmain.f90 revision 57 by guez, Mon Jan 30 12:54:02 2012 UTC trunk/phylmd/clmain.f revision 82 by guez, Wed Mar 5 14:57:53 2014 UTC
# Line 4  module clmain_m Line 4  module clmain_m
4    
5  contains  contains
6    
7    SUBROUTINE clmain(dtime, itap, date0, pctsrf, pctsrf_new, t, q, u, v,&    SUBROUTINE clmain(dtime, itap, pctsrf, pctsrf_new, t, q, u, v, &
8         jour, rmu0, co2_ppm, ok_veget, ocean, npas, nexca, ts,&         jour, rmu0, co2_ppm, ok_veget, ocean, ts, &
9         soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil,&         soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, &
10         qsol, paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat,&         qsol, paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat, &
11         rain_f, snow_f, solsw, sollw, sollwdown, fder, rlon, rlat, cufi,&         rain_fall, snow_f, solsw, sollw, fder, rlon, rlat, &
12         cvfi, rugos, debut, lafin, agesno, rugoro, d_t, d_q, d_u, d_v,&         rugos, debut, agesno, rugoro, d_t, d_q, d_u, d_v, &
13         d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2,&         d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, &
14         dflux_t, dflux_q, zcoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh,&         dflux_t, dflux_q, ycoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh, &
15         capcl, oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl,&         capcl, oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl, &
16         fqcalving, ffonte, run_off_lic_0, flux_o, flux_g, tslab, seaice)         fqcalving, ffonte, run_off_lic_0, flux_o, flux_g, tslab, seaice)
17    
18      ! From phylmd/clmain.F, version 1.6 2005/11/16 14:47:19      ! From phylmd/clmain.F, version 1.6 2005/11/16 14:47:19
19      ! Author: Z.X. Li (LMD/CNRS), date: 1993/08/18      ! Author: Z. X. Li (LMD/CNRS), date: 1993/08/18
20      ! Objet : interface de "couche limite" (diffusion verticale)      ! Objet : interface de couche limite (diffusion verticale)
21    
22      ! Tout ce qui a trait aux traceurs est dans "phytrac" maintenant.      ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
23      ! Pour l'instant le calcul de la couche limite pour les traceurs      ! de la couche limite pour les traceurs se fait avec "cltrac" et
24      ! se fait avec "cltrac" et ne tient pas compte de la différentiation      ! ne tient pas compte de la différentiation des sous-fractions de
25      ! des sous-fractions de sol.      ! sol.
26    
27      ! Pour pouvoir extraire les coefficients d'échanges et le vent      ! Pour pouvoir extraire les coefficients d'échanges et le vent
28      ! dans la première couche, trois champs supplémentaires ont été      ! dans la première couche, trois champs ont été créés : "ycoefh",
29      ! créés : "zcoefh", "zu1" et "zv1". Pour l'instant nous avons      ! "zu1" et "zv1". Nous avons moyenné les valeurs de ces trois
30      ! moyenné les valeurs de ces trois champs sur les 4 sous-surfaces      ! champs sur les quatre sous-surfaces du modèle.
     ! du modèle. Dans l'avenir, si les informations des sous-surfaces  
     ! doivent être prises en compte, il faudra sortir ces mêmes champs  
     ! en leur ajoutant une dimension, c'est-à-dire "nbsrf" (nombre de  
     ! sous-surfaces).  
31    
32      use calendar, ONLY : ymds2ju      use calendar, ONLY: ymds2ju
33      use clqh_m, only: clqh      use clqh_m, only: clqh
34        use clvent_m, only: clvent
35      use coefkz_m, only: coefkz      use coefkz_m, only: coefkz
36      use coefkzmin_m, only: coefkzmin      use coefkzmin_m, only: coefkzmin
37      USE conf_phys_m, ONLY : iflag_pbl      USE conf_gcm_m, ONLY: prt_level
38      USE dimens_m, ONLY : iim, jjm      USE conf_phys_m, ONLY: iflag_pbl
39      USE dimphy, ONLY : klev, klon, zmasq      USE dimens_m, ONLY: iim, jjm
40      USE dimsoil, ONLY : nsoilmx      USE dimphy, ONLY: klev, klon, zmasq
41      USE dynetat0_m, ONLY : day_ini      USE dimsoil, ONLY: nsoilmx
42      USE gath_cpl, ONLY : gath2cpl      USE dynetat0_m, ONLY: day_ini
43        USE gath_cpl, ONLY: gath2cpl
44      use hbtm_m, only: hbtm      use hbtm_m, only: hbtm
45      USE histcom, ONLY : histbeg_totreg, histdef, histend, histsync      USE histbeg_totreg_m, ONLY: histbeg_totreg
46        USE histdef_m, ONLY: histdef
47        USE histend_m, ONLY: histend
48        USE histsync_m, ONLY: histsync
49      use histwrite_m, only: histwrite      use histwrite_m, only: histwrite
50      USE indicesol, ONLY : epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf      USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
51      USE conf_gcm_m, ONLY : prt_level      USE suphec_m, ONLY: rd, rg, rkappa
52      USE suphec_m, ONLY : rd, rg, rkappa      USE temps, ONLY: annee_ref, itau_phy
53      USE temps, ONLY : annee_ref, itau_phy      use ustarhb_m, only: ustarhb
54        use vdif_kcay_m, only: vdif_kcay
55      use yamada4_m, only: yamada4      use yamada4_m, only: yamada4
56    
57      ! Arguments:      ! Arguments:
58    
59      REAL, INTENT (IN) :: dtime ! interval du temps (secondes)      REAL, INTENT(IN):: dtime ! interval du temps (secondes)
60      REAL date0      INTEGER, INTENT(IN):: itap ! numero du pas de temps
61      ! date0----input-R- jour initial      REAL, INTENT(inout):: pctsrf(klon, nbsrf)
62      INTEGER, INTENT (IN) :: itap  
63      ! itap-----input-I- numero du pas de temps      ! la nouvelle repartition des surfaces sortie de l'interface
64      REAL, INTENT(IN):: t(klon, klev), q(klon, klev)      REAL, INTENT(out):: pctsrf_new(klon, nbsrf)
65      ! t--------input-R- temperature (K)  
66      ! q--------input-R- vapeur d'eau (kg/kg)      REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
67      REAL, INTENT (IN):: u(klon, klev), v(klon, klev)      REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg/kg)
68      ! u--------input-R- vitesse u      REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
69      ! v--------input-R- vitesse v      INTEGER, INTENT(IN):: jour ! jour de l'annee en cours
70      REAL, INTENT (IN):: paprs(klon, klev+1)      REAL, intent(in):: rmu0(klon) ! cosinus de l'angle solaire zenithal    
71      ! paprs----input-R- pression a intercouche (Pa)      REAL co2_ppm ! taux CO2 atmosphere
72      REAL, INTENT (IN):: pplay(klon, klev)      LOGICAL ok_veget
73      ! pplay----input-R- pression au milieu de couche (Pa)      CHARACTER(len=*), INTENT(IN):: ocean
74      REAL, INTENT (IN):: rlon(klon), rlat(klon)      REAL ts(klon, nbsrf) ! input-R- temperature du sol (en Kelvin)
75      ! rlat-----input-R- latitude en degree      LOGICAL, INTENT(IN):: soil_model
76      REAL cufi(klon), cvfi(klon)      REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
77      ! cufi-----input-R- resolution des mailles en x (m)      REAL ksta, ksta_ter
78      ! cvfi-----input-R- resolution des mailles en y (m)      LOGICAL ok_kzmin
79        REAL ftsoil(klon, nsoilmx, nbsrf)
80        REAL qsol(klon)
81        REAL, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)
82        REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
83        REAL snow(klon, nbsrf)
84        REAL qsurf(klon, nbsrf)
85        REAL evap(klon, nbsrf)
86        REAL albe(klon, nbsrf)
87        REAL alblw(klon, nbsrf)
88    
89        REAL fluxlat(klon, nbsrf)
90    
91        REAL, intent(in):: rain_fall(klon), snow_f(klon)
92        REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)
93        REAL fder(klon)
94        REAL, INTENT(IN):: rlon(klon)
95        REAL, INTENT(IN):: rlat(klon) ! latitude en degrés
96    
97        REAL rugos(klon, nbsrf)
98        ! rugos----input-R- longeur de rugosite (en m)
99    
100        LOGICAL, INTENT(IN):: debut
101        real agesno(klon, nbsrf)
102        REAL, INTENT(IN):: rugoro(klon)
103    
104      REAL d_t(klon, klev), d_q(klon, klev)      REAL d_t(klon, klev), d_q(klon, klev)
105      ! d_t------output-R- le changement pour "t"      ! d_t------output-R- le changement pour "t"
106      ! d_q------output-R- le changement pour "q"      ! d_q------output-R- le changement pour "q"
107      REAL d_u(klon, klev), d_v(klon, klev)  
108      ! d_u------output-R- le changement pour "u"      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
109      ! d_v------output-R- le changement pour "v"      ! changement pour "u" et "v"
110    
111        REAL d_ts(klon, nbsrf)
112        ! d_ts-----output-R- le changement pour "ts"
113    
114      REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)      REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)
115      ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)      ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
116      !                    (orientation positive vers le bas)      !                    (orientation positive vers le bas)
117      ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)      ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
118    
119        REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)
120        ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
121        ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
122    
123        REAL, INTENT(out):: cdragh(klon), cdragm(klon)
124        real q2(klon, klev+1, nbsrf)
125    
126      REAL dflux_t(klon), dflux_q(klon)      REAL dflux_t(klon), dflux_q(klon)
127      ! dflux_t derive du flux sensible      ! dflux_t derive du flux sensible
128      ! dflux_q derive du flux latent      ! dflux_q derive du flux latent
129      !IM "slab" ocean      !IM "slab" ocean
130      REAL flux_o(klon), flux_g(klon)  
131      !IM "slab" ocean      REAL, intent(out):: ycoefh(klon, klev)
132      ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')      REAL, intent(out):: zu1(klon)
133      ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')      REAL zv1(klon)
134      REAL y_flux_o(klon), y_flux_g(klon)      REAL t2m(klon, nbsrf), q2m(klon, nbsrf)
135      REAL tslab(klon), ytslab(klon)      REAL u10m(klon, nbsrf), v10m(klon, nbsrf)
136      ! tslab-in/output-R temperature du slab ocean (en Kelvin)  
137      ! uniqmnt pour slab      !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds
138      REAL seaice(klon), y_seaice(klon)      ! physiq ce qui permet de sortir les grdeurs par sous surface)
139      ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')      REAL pblh(klon, nbsrf)
140      REAL y_fqcalving(klon), y_ffonte(klon)      ! pblh------- HCL
141        REAL capcl(klon, nbsrf)
142        REAL oliqcl(klon, nbsrf)
143        REAL cteicl(klon, nbsrf)
144        REAL pblt(klon, nbsrf)
145        ! pblT------- T au nveau HCL
146        REAL therm(klon, nbsrf)
147        REAL trmb1(klon, nbsrf)
148        ! trmb1-------deep_cape
149        REAL trmb2(klon, nbsrf)
150        ! trmb2--------inhibition
151        REAL trmb3(klon, nbsrf)
152        ! trmb3-------Point Omega
153        REAL plcl(klon, nbsrf)
154      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
155      ! ffonte----Flux thermique utilise pour fondre la neige      ! ffonte----Flux thermique utilise pour fondre la neige
156      ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la      ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
157      !           hauteur de neige, en kg/m2/s      !           hauteur de neige, en kg/m2/s
158      REAL run_off_lic_0(klon), y_run_off_lic_0(klon)      REAL run_off_lic_0(klon)
   
     REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)  
     ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal  
     ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal  
     REAL rugmer(klon), agesno(klon, nbsrf)  
     REAL, INTENT (IN) :: rugoro(klon)  
     REAL cdragh(klon), cdragm(klon)  
     ! jour de l'annee en cours                  
     INTEGER jour  
     REAL rmu0(klon) ! cosinus de l'angle solaire zenithal      
     ! taux CO2 atmosphere                      
     REAL co2_ppm  
     LOGICAL, INTENT (IN) :: debut  
     LOGICAL, INTENT (IN) :: lafin  
     LOGICAL ok_veget  
     CHARACTER (len=*), INTENT (IN) :: ocean  
     INTEGER npas, nexca  
   
     REAL pctsrf(klon, nbsrf)  
     REAL ts(klon, nbsrf)  
     ! ts-------input-R- temperature du sol (en Kelvin)  
     REAL d_ts(klon, nbsrf)  
     ! d_ts-----output-R- le changement pour "ts"  
     REAL snow(klon, nbsrf)  
     REAL qsurf(klon, nbsrf)  
     REAL evap(klon, nbsrf)  
     REAL albe(klon, nbsrf)  
     REAL alblw(klon, nbsrf)  
159    
160      REAL fluxlat(klon, nbsrf)      REAL flux_o(klon), flux_g(klon)
161        !IM "slab" ocean
162        ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')
163        ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')
164    
165      REAL rain_f(klon), snow_f(klon)      REAL tslab(klon)
166      REAL fder(klon)      ! tslab-in/output-R temperature du slab ocean (en Kelvin)
167        ! uniqmnt pour slab
168    
169      REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)      REAL seaice(klon)
170      REAL rugos(klon, nbsrf)      ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')
     ! rugos----input-R- longeur de rugosite (en m)  
     ! la nouvelle repartition des surfaces sortie de l'interface  
     REAL pctsrf_new(klon, nbsrf)  
171    
172      REAL zcoefh(klon, klev)      ! Local:
     REAL zu1(klon)  
     REAL zv1(klon)  
173    
174      !$$$ PB ajout pour soil      REAL y_flux_o(klon), y_flux_g(klon)
175      LOGICAL, INTENT (IN) :: soil_model      real ytslab(klon)
176      !IM ajout seuils cdrm, cdrh      real y_seaice(klon)
177      REAL cdmmax, cdhmax      REAL y_fqcalving(klon), y_ffonte(klon)
178        real y_run_off_lic_0(klon)
179    
180      REAL ksta, ksta_ter      REAL rugmer(klon)
     LOGICAL ok_kzmin  
181    
     REAL ftsoil(klon, nsoilmx, nbsrf)  
182      REAL ytsoil(klon, nsoilmx)      REAL ytsoil(klon, nsoilmx)
     REAL qsol(klon)  
   
     EXTERNAL clvent, calbeta, cltrac  
183    
184      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
185      REAL yalb(klon)      REAL yalb(klon)
# Line 170  contains Line 189  contains
189      ! la premiere couche      ! la premiere couche
190      REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)      REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)
191      REAL yrain_f(klon), ysnow_f(klon)      REAL yrain_f(klon), ysnow_f(klon)
192      REAL ysollw(klon), ysolsw(klon), ysollwdown(klon)      REAL ysollw(klon), ysolsw(klon)
193      REAL yfder(klon), ytaux(klon), ytauy(klon)      REAL yfder(klon), ytaux(klon), ytauy(klon)
194      REAL yrugm(klon), yrads(klon), yrugoro(klon)      REAL yrugm(klon), yrads(klon), yrugoro(klon)
195    
# Line 182  contains Line 201  contains
201      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)
202      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)
203      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
204      REAL ycoefh(klon, klev), ycoefm(klon, klev)      REAL coefh(klon, klev), coefm(klon, klev)
205      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
206      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
207      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)
208    
     LOGICAL ok_nonloc  
     PARAMETER (ok_nonloc=.FALSE.)  
209      REAL ycoefm0(klon, klev), ycoefh0(klon, klev)      REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
210    
211      REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)      REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
212      REAL ykmm(klon, klev+1), ykmn(klon, klev+1)      REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
213      REAL ykmq(klon, klev+1)      REAL ykmq(klon, klev+1)
214      REAL yq2(klon, klev+1), q2(klon, klev+1, nbsrf)      REAL yq2(klon, klev+1)
215      REAL q2diag(klon, klev+1)      REAL q2diag(klon, klev+1)
216    
217      REAL u1lay(klon), v1lay(klon)      REAL u1lay(klon), v1lay(klon)
# Line 211  contains Line 228  contains
228    
229      ! maf pour sorties IOISPL en cas de debugagage      ! maf pour sorties IOISPL en cas de debugagage
230    
231      CHARACTER (80) cldebug      CHARACTER(80) cldebug
232      SAVE cldebug      SAVE cldebug
233      CHARACTER (8) cl_surf(nbsrf)      CHARACTER(8) cl_surf(nbsrf)
234      SAVE cl_surf      SAVE cl_surf
235      INTEGER nhoridbg, nidbg      INTEGER nhoridbg, nidbg
236      SAVE nhoridbg, nidbg      SAVE nhoridbg, nidbg
# Line 224  contains Line 241  contains
241      LOGICAL first_appel      LOGICAL first_appel
242      SAVE first_appel      SAVE first_appel
243      DATA first_appel/ .TRUE./      DATA first_appel/ .TRUE./
244      LOGICAL :: debugindex = .FALSE.      LOGICAL:: debugindex = .FALSE.
245      INTEGER idayref      INTEGER idayref
     REAL t2m(klon, nbsrf), q2m(klon, nbsrf)  
     REAL u10m(klon, nbsrf), v10m(klon, nbsrf)  
246    
247      REAL yt2m(klon), yq2m(klon), yu10m(klon)      REAL yt2m(klon), yq2m(klon), yu10m(klon)
248      REAL yustar(klon)      REAL yustar(klon)
# Line 238  contains Line 253  contains
253      ! -- LOOP      ! -- LOOP
254    
255      REAL yt10m(klon), yq10m(klon)      REAL yt10m(klon), yq10m(klon)
     !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds  
     ! physiq ce qui permet de sortir les grdeurs par sous surface)  
     REAL pblh(klon, nbsrf)  
     ! pblh------- HCL  
     REAL plcl(klon, nbsrf)  
     REAL capcl(klon, nbsrf)  
     REAL oliqcl(klon, nbsrf)  
     REAL cteicl(klon, nbsrf)  
     REAL pblt(klon, nbsrf)  
     ! pblT------- T au nveau HCL  
     REAL therm(klon, nbsrf)  
     REAL trmb1(klon, nbsrf)  
     ! trmb1-------deep_cape  
     REAL trmb2(klon, nbsrf)  
     ! trmb2--------inhibition  
     REAL trmb3(klon, nbsrf)  
     ! trmb3-------Point Omega  
256      REAL ypblh(klon)      REAL ypblh(klon)
257      REAL ylcl(klon)      REAL ylcl(klon)
258      REAL ycapcl(klon)      REAL ycapcl(klon)
# Line 265  contains Line 263  contains
263      REAL ytrmb1(klon)      REAL ytrmb1(klon)
264      REAL ytrmb2(klon)      REAL ytrmb2(klon)
265      REAL ytrmb3(klon)      REAL ytrmb3(klon)
     REAL y_cd_h(klon), y_cd_m(klon)  
266      REAL uzon(klon), vmer(klon)      REAL uzon(klon), vmer(klon)
267      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
268      REAL psfce(klon), patm(klon)      REAL psfce(klon), patm(klon)
# Line 277  contains Line 274  contains
274      LOGICAL zxli      LOGICAL zxli
275      PARAMETER (zxli=.FALSE.)      PARAMETER (zxli=.FALSE.)
276    
     REAL zt, zqs, zdelta, zcor  
     REAL t_coup  
     PARAMETER (t_coup=273.15)  
   
     CHARACTER (len=20) :: modname = 'clmain'  
   
277      !------------------------------------------------------------      !------------------------------------------------------------
278    
279      ytherm = 0.      ytherm = 0.
# Line 349  contains Line 340  contains
340      ytauy = 0.      ytauy = 0.
341      ysolsw = 0.      ysolsw = 0.
342      ysollw = 0.      ysollw = 0.
     ysollwdown = 0.  
343      yrugos = 0.      yrugos = 0.
344      yu1 = 0.      yu1 = 0.
345      yv1 = 0.      yv1 = 0.
# Line 385  contains Line 375  contains
375      d_q = 0.      d_q = 0.
376      d_u = 0.      d_u = 0.
377      d_v = 0.      d_v = 0.
378      zcoefh = 0.      ycoefh = 0.
379    
380      ! Boucler sur toutes les sous-fractions du sol:      ! Boucler sur toutes les sous-fractions du sol:
381    
# Line 422  contains Line 412  contains
412            CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)            CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)
413         END IF         END IF
414    
415         IF (knon == 0) CYCLE         if_knon: IF (knon /= 0) then
   
        DO j = 1, knon  
           i = ni(j)  
           ypct(j) = pctsrf(i, nsrf)  
           yts(j) = ts(i, nsrf)  
           ytslab(i) = tslab(i)  
           ysnow(j) = snow(i, nsrf)  
           yqsurf(j) = qsurf(i, nsrf)  
           yalb(j) = albe(i, nsrf)  
           yalblw(j) = alblw(i, nsrf)  
           yrain_f(j) = rain_f(i)  
           ysnow_f(j) = snow_f(i)  
           yagesno(j) = agesno(i, nsrf)  
           yfder(j) = fder(i)  
           ytaux(j) = flux_u(i, 1, nsrf)  
           ytauy(j) = flux_v(i, 1, nsrf)  
           ysolsw(j) = solsw(i, nsrf)  
           ysollw(j) = sollw(i, nsrf)  
           ysollwdown(j) = sollwdown(i)  
           yrugos(j) = rugos(i, nsrf)  
           yrugoro(j) = rugoro(i)  
           yu1(j) = u1lay(i)  
           yv1(j) = v1lay(i)  
           yrads(j) = ysolsw(j) + ysollw(j)  
           ypaprs(j, klev+1) = paprs(i, klev+1)  
           y_run_off_lic_0(j) = run_off_lic_0(i)  
           yu10mx(j) = u10m(i, nsrf)  
           yu10my(j) = v10m(i, nsrf)  
           ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))  
        END DO  
   
        ! IF bucket model for continent, copy soil water content  
        IF (nsrf == is_ter .AND. .NOT. ok_veget) THEN  
           DO j = 1, knon  
              i = ni(j)  
              yqsol(j) = qsol(i)  
           END DO  
        ELSE  
           yqsol = 0.  
        END IF  
        !$$$ PB ajour pour soil  
        DO k = 1, nsoilmx  
416            DO j = 1, knon            DO j = 1, knon
417               i = ni(j)               i = ni(j)
418               ytsoil(j, k) = ftsoil(i, k, nsrf)               ypct(j) = pctsrf(i, nsrf)
419                 yts(j) = ts(i, nsrf)
420                 ytslab(i) = tslab(i)
421                 ysnow(j) = snow(i, nsrf)
422                 yqsurf(j) = qsurf(i, nsrf)
423                 yalb(j) = albe(i, nsrf)
424                 yalblw(j) = alblw(i, nsrf)
425                 yrain_f(j) = rain_fall(i)
426                 ysnow_f(j) = snow_f(i)
427                 yagesno(j) = agesno(i, nsrf)
428                 yfder(j) = fder(i)
429                 ytaux(j) = flux_u(i, 1, nsrf)
430                 ytauy(j) = flux_v(i, 1, nsrf)
431                 ysolsw(j) = solsw(i, nsrf)
432                 ysollw(j) = sollw(i, nsrf)
433                 yrugos(j) = rugos(i, nsrf)
434                 yrugoro(j) = rugoro(i)
435                 yu1(j) = u1lay(i)
436                 yv1(j) = v1lay(i)
437                 yrads(j) = ysolsw(j) + ysollw(j)
438                 ypaprs(j, klev+1) = paprs(i, klev+1)
439                 y_run_off_lic_0(j) = run_off_lic_0(i)
440                 yu10mx(j) = u10m(i, nsrf)
441                 yu10my(j) = v10m(i, nsrf)
442                 ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))
443            END DO            END DO
        END DO  
        DO k = 1, klev  
           DO j = 1, knon  
              i = ni(j)  
              ypaprs(j, k) = paprs(i, k)  
              ypplay(j, k) = pplay(i, k)  
              ydelp(j, k) = delp(i, k)  
              yu(j, k) = u(i, k)  
              yv(j, k) = v(i, k)  
              yt(j, k) = t(i, k)  
              yq(j, k) = q(i, k)  
           END DO  
        END DO  
444    
445         ! calculer Cdrag et les coefficients d'echange            ! IF bucket model for continent, copy soil water content
446         CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,&            IF (nsrf == is_ter .AND. .NOT. ok_veget) THEN
447              yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)               DO j = 1, knon
448         IF (iflag_pbl == 1) THEN                  i = ni(j)
449            CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)                  yqsol(j) = qsol(i)
           DO k = 1, klev  
              DO i = 1, knon  
                 ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))  
                 ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))  
450               END DO               END DO
451            END DO            ELSE
452         END IF               yqsol = 0.
453              END IF
        ! on seuille ycoefm et ycoefh  
        IF (nsrf == is_oce) THEN  
           DO j = 1, knon  
              ycoefm(j, 1) = min(ycoefm(j, 1), cdmmax)  
              ycoefh(j, 1) = min(ycoefh(j, 1), cdhmax)  
           END DO  
        END IF  
   
        IF (ok_kzmin) THEN  
           ! Calcul d'une diffusion minimale pour les conditions tres stables  
           CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, ycoefm(:, 1), &  
                ycoefm0, ycoefh0)  
454    
455            DO k = 1, klev            DO k = 1, nsoilmx
456               DO i = 1, knon               DO j = 1, knon
457                  ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))                  i = ni(j)
458                  ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))                  ytsoil(j, k) = ftsoil(i, k, nsrf)
459               END DO               END DO
460            END DO            END DO
        END IF  
461    
        IF (iflag_pbl >= 3) THEN  
           ! MELLOR ET YAMADA adapté à Mars, Richard Fournier et Frédéric Hourdin  
           yzlay(1:knon, 1) = rd*yt(1:knon, 1)/(0.5*(ypaprs(1:knon, &  
                1)+ypplay(1:knon, 1)))*(ypaprs(1:knon, 1)-ypplay(1:knon, 1))/rg  
           DO k = 2, klev  
              yzlay(1:knon, k) = yzlay(1:knon, k-1) &  
                   + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &  
                   / ypaprs(1:knon, k) &  
                   * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg  
           END DO  
462            DO k = 1, klev            DO k = 1, klev
              yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &  
                   / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))  
           END DO  
           yzlev(1:knon, 1) = 0.  
           yzlev(1:knon, klev+1) = 2.*yzlay(1:knon, klev) - yzlay(1:knon, klev-1)  
           DO k = 2, klev  
              yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))  
           END DO  
           DO k = 1, klev + 1  
463               DO j = 1, knon               DO j = 1, knon
464                  i = ni(j)                  i = ni(j)
465                  yq2(j, k) = q2(i, k, nsrf)                  ypaprs(j, k) = paprs(i, k)
466                    ypplay(j, k) = pplay(i, k)
467                    ydelp(j, k) = delp(i, k)
468                    yu(j, k) = u(i, k)
469                    yv(j, k) = v(i, k)
470                    yt(j, k) = t(i, k)
471                    yq(j, k) = q(i, k)
472               END DO               END DO
473            END DO            END DO
474    
475            y_cd_m(1:knon) = ycoefm(1:knon, 1)            ! calculer Cdrag et les coefficients d'echange
476            y_cd_h(1:knon) = ycoefh(1:knon, 1)            CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &
477            CALL ustarhb(knon, yu, yv, y_cd_m, yustar)                 yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))
478              IF (iflag_pbl == 1) THEN
479                 CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
480                 coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
481                 coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
482              END IF
483    
484            IF (prt_level>9) THEN            ! on met un seuil pour coefm et coefh
485               PRINT *, 'USTAR = ', yustar            IF (nsrf == is_oce) THEN
486                 coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
487                 coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
488            END IF            END IF
489    
490            ! iflag_pbl peut être utilisé comme longueur de mélange            IF (ok_kzmin) THEN
491                 ! Calcul d'une diffusion minimale pour les conditions tres stables
492                 CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
493                      coefm(:knon, 1), ycoefm0, ycoefh0)
494                 coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
495                 coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
496               END IF
497    
498              IF (iflag_pbl >= 3) THEN
499                 ! MELLOR ET YAMADA adapté à Mars, Richard Fournier et
500                 ! Frédéric Hourdin
501                 yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
502                      + ypplay(:knon, 1))) &
503                      * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
504                 DO k = 2, klev
505                    yzlay(1:knon, k) = yzlay(1:knon, k-1) &
506                         + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
507                         / ypaprs(1:knon, k) &
508                         * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
509                 END DO
510                 DO k = 1, klev
511                    yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
512                         / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
513                 END DO
514                 yzlev(1:knon, 1) = 0.
515                 yzlev(:knon, klev+1) = 2. * yzlay(:knon, klev) &
516                      - yzlay(:knon, klev - 1)
517                 DO k = 2, klev
518                    yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
519                 END DO
520                 DO k = 1, klev + 1
521                    DO j = 1, knon
522                       i = ni(j)
523                       yq2(j, k) = q2(i, k, nsrf)
524                    END DO
525                 END DO
526    
527                 CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
528    
529            IF (iflag_pbl >= 11) THEN               IF (prt_level > 9) THEN
530               CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &                  PRINT *, 'USTAR = ', yustar
531                    yu, yv, yteta, y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, &               END IF
532                    iflag_pbl)  
533            ELSE               ! iflag_pbl peut être utilisé comme longueur de mélange
534               CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &  
535                    y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)               IF (iflag_pbl >= 11) THEN
536                    CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &
537                         yu, yv, yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, &
538                         yustar, iflag_pbl)
539                 ELSE
540                    CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
541                         coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
542                 END IF
543    
544                 coefm(:knon, 2:) = ykmm(:knon, 2:klev)
545                 coefh(:knon, 2:) = ykmn(:knon, 2:klev)
546            END IF            END IF
547    
548            ycoefm(1:knon, 1) = y_cd_m(1:knon)            ! calculer la diffusion des vitesses "u" et "v"
549            ycoefh(1:knon, 1) = y_cd_h(1:knon)            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &
550            ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)                 ypplay, ydelp, y_d_u, y_flux_u)
551            ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
552         END IF                 ypplay, ydelp, y_d_v, y_flux_v)
553    
554              ! pour le couplage
555              ytaux = y_flux_u(:, 1)
556              ytauy = y_flux_v(:, 1)
557    
558              ! calculer la diffusion de "q" et de "h"
559              CALL clqh(dtime, itap, jour, debut, rlat, knon, nsrf, ni, pctsrf, &
560                   soil_model, ytsoil, yqsol, ok_veget, ocean, rmu0, co2_ppm, &
561                   yrugos, yrugoro, yu1, yv1, coefh(:knon, :), yt, yq, yts, &
562                   ypaprs, ypplay, ydelp, yrads, yalb, yalblw, ysnow, yqsurf, &
563                   yrain_f, ysnow_f, yfder, ysolsw, yfluxlat, pctsrf_new, &
564                   yagesno, y_d_t, y_d_q, y_d_ts, yz0_new, y_flux_t, y_flux_q, &
565                   y_dflux_t, y_dflux_q, y_fqcalving, y_ffonte, y_run_off_lic_0, &
566                   y_flux_o, y_flux_g, ytslab, y_seaice)
567    
568         ! calculer la diffusion des vitesses "u" et "v"            ! calculer la longueur de rugosite sur ocean
569         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &            yrugm = 0.
570              ydelp, y_d_u, y_flux_u)            IF (nsrf == is_oce) THEN
571         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &               DO j = 1, knon
572              ydelp, y_d_v, y_flux_v)                  yrugm(j) = 0.018*coefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
573                         0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))
574         ! pour le couplage                  yrugm(j) = max(1.5E-05, yrugm(j))
575         ytaux = y_flux_u(:, 1)               END DO
576         ytauy = y_flux_v(:, 1)            END IF
   
        ! calculer la diffusion de "q" et de "h"  
        CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,&  
             cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil,&  
             yqsol, ok_veget, ocean, npas, nexca, rmu0, co2_ppm, yrugos,&  
             yrugoro, yu1, yv1, ycoefh, yt, yq, yts, ypaprs, ypplay,&  
             ydelp, yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, &  
             yfder, ytaux, ytauy, ywindsp, ysollw, ysollwdown, ysolsw,&  
             yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts,&  
             yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,&  
             y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g,&  
             ytslab, y_seaice)  
   
        ! calculer la longueur de rugosite sur ocean  
        yrugm = 0.  
        IF (nsrf == is_oce) THEN  
577            DO j = 1, knon            DO j = 1, knon
578               yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &               y_dflux_t(j) = y_dflux_t(j)*ypct(j)
579                    0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2))               y_dflux_q(j) = y_dflux_q(j)*ypct(j)
580               yrugm(j) = max(1.5E-05, yrugm(j))               yu1(j) = yu1(j)*ypct(j)
581                 yv1(j) = yv1(j)*ypct(j)
582            END DO            END DO
        END IF  
        DO j = 1, knon  
           y_dflux_t(j) = y_dflux_t(j)*ypct(j)  
           y_dflux_q(j) = y_dflux_q(j)*ypct(j)  
           yu1(j) = yu1(j)*ypct(j)  
           yv1(j) = yv1(j)*ypct(j)  
        END DO  
583    
584         DO k = 1, klev            DO k = 1, klev
585            DO j = 1, knon               DO j = 1, knon
586               i = ni(j)                  i = ni(j)
587               ycoefh(j, k) = ycoefh(j, k)*ypct(j)                  coefh(j, k) = coefh(j, k)*ypct(j)
588               ycoefm(j, k) = ycoefm(j, k)*ypct(j)                  coefm(j, k) = coefm(j, k)*ypct(j)
589               y_d_t(j, k) = y_d_t(j, k)*ypct(j)                  y_d_t(j, k) = y_d_t(j, k)*ypct(j)
590               y_d_q(j, k) = y_d_q(j, k)*ypct(j)                  y_d_q(j, k) = y_d_q(j, k)*ypct(j)
591               flux_t(i, k, nsrf) = y_flux_t(j, k)                  flux_t(i, k, nsrf) = y_flux_t(j, k)
592               flux_q(i, k, nsrf) = y_flux_q(j, k)                  flux_q(i, k, nsrf) = y_flux_q(j, k)
593               flux_u(i, k, nsrf) = y_flux_u(j, k)                  flux_u(i, k, nsrf) = y_flux_u(j, k)
594               flux_v(i, k, nsrf) = y_flux_v(j, k)                  flux_v(i, k, nsrf) = y_flux_v(j, k)
595               y_d_u(j, k) = y_d_u(j, k)*ypct(j)                  y_d_u(j, k) = y_d_u(j, k)*ypct(j)
596               y_d_v(j, k) = y_d_v(j, k)*ypct(j)                  y_d_v(j, k) = y_d_v(j, k)*ypct(j)
597                 END DO
598            END DO            END DO
        END DO  
599    
600         evap(:, nsrf) = -flux_q(:, 1, nsrf)            evap(:, nsrf) = -flux_q(:, 1, nsrf)
601    
602         albe(:, nsrf) = 0.            albe(:, nsrf) = 0.
603         alblw(:, nsrf) = 0.            alblw(:, nsrf) = 0.
604         snow(:, nsrf) = 0.            snow(:, nsrf) = 0.
605         qsurf(:, nsrf) = 0.            qsurf(:, nsrf) = 0.
606         rugos(:, nsrf) = 0.            rugos(:, nsrf) = 0.
607         fluxlat(:, nsrf) = 0.            fluxlat(:, nsrf) = 0.
        DO j = 1, knon  
           i = ni(j)  
           d_ts(i, nsrf) = y_d_ts(j)  
           albe(i, nsrf) = yalb(j)  
           alblw(i, nsrf) = yalblw(j)  
           snow(i, nsrf) = ysnow(j)  
           qsurf(i, nsrf) = yqsurf(j)  
           rugos(i, nsrf) = yz0_new(j)  
           fluxlat(i, nsrf) = yfluxlat(j)  
           IF (nsrf == is_oce) THEN  
              rugmer(i) = yrugm(j)  
              rugos(i, nsrf) = yrugm(j)  
           END IF  
           agesno(i, nsrf) = yagesno(j)  
           fqcalving(i, nsrf) = y_fqcalving(j)  
           ffonte(i, nsrf) = y_ffonte(j)  
           cdragh(i) = cdragh(i) + ycoefh(j, 1)  
           cdragm(i) = cdragm(i) + ycoefm(j, 1)  
           dflux_t(i) = dflux_t(i) + y_dflux_t(j)  
           dflux_q(i) = dflux_q(i) + y_dflux_q(j)  
           zu1(i) = zu1(i) + yu1(j)  
           zv1(i) = zv1(i) + yv1(j)  
        END DO  
        IF (nsrf == is_ter) THEN  
608            DO j = 1, knon            DO j = 1, knon
609               i = ni(j)               i = ni(j)
610               qsol(i) = yqsol(j)               d_ts(i, nsrf) = y_d_ts(j)
611                 albe(i, nsrf) = yalb(j)
612                 alblw(i, nsrf) = yalblw(j)
613                 snow(i, nsrf) = ysnow(j)
614                 qsurf(i, nsrf) = yqsurf(j)
615                 rugos(i, nsrf) = yz0_new(j)
616                 fluxlat(i, nsrf) = yfluxlat(j)
617                 IF (nsrf == is_oce) THEN
618                    rugmer(i) = yrugm(j)
619                    rugos(i, nsrf) = yrugm(j)
620                 END IF
621                 agesno(i, nsrf) = yagesno(j)
622                 fqcalving(i, nsrf) = y_fqcalving(j)
623                 ffonte(i, nsrf) = y_ffonte(j)
624                 cdragh(i) = cdragh(i) + coefh(j, 1)
625                 cdragm(i) = cdragm(i) + coefm(j, 1)
626                 dflux_t(i) = dflux_t(i) + y_dflux_t(j)
627                 dflux_q(i) = dflux_q(i) + y_dflux_q(j)
628                 zu1(i) = zu1(i) + yu1(j)
629                 zv1(i) = zv1(i) + yv1(j)
630            END DO            END DO
631         END IF            IF (nsrf == is_ter) THEN
632         IF (nsrf == is_lic) THEN               DO j = 1, knon
633                    i = ni(j)
634                    qsol(i) = yqsol(j)
635                 END DO
636              END IF
637              IF (nsrf == is_lic) THEN
638                 DO j = 1, knon
639                    i = ni(j)
640                    run_off_lic_0(i) = y_run_off_lic_0(j)
641                 END DO
642              END IF
643              !$$$ PB ajout pour soil
644              ftsoil(:, :, nsrf) = 0.
645              DO k = 1, nsoilmx
646                 DO j = 1, knon
647                    i = ni(j)
648                    ftsoil(i, k, nsrf) = ytsoil(j, k)
649                 END DO
650              END DO
651    
652            DO j = 1, knon            DO j = 1, knon
653               i = ni(j)               i = ni(j)
654               run_off_lic_0(i) = y_run_off_lic_0(j)               DO k = 1, klev
655                    d_t(i, k) = d_t(i, k) + y_d_t(j, k)
656                    d_q(i, k) = d_q(i, k) + y_d_q(j, k)
657                    d_u(i, k) = d_u(i, k) + y_d_u(j, k)
658                    d_v(i, k) = d_v(i, k) + y_d_v(j, k)
659                    ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
660                 END DO
661            END DO            END DO
662         END IF  
663         !$$$ PB ajout pour soil            !cc diagnostic t, q a 2m et u, v a 10m
664         ftsoil(:, :, nsrf) = 0.  
        DO k = 1, nsoilmx  
665            DO j = 1, knon            DO j = 1, knon
666               i = ni(j)               i = ni(j)
667               ftsoil(i, k, nsrf) = ytsoil(j, k)               uzon(j) = yu(j, 1) + y_d_u(j, 1)
668            END DO               vmer(j) = yv(j, 1) + y_d_v(j, 1)
669         END DO               tair1(j) = yt(j, 1) + y_d_t(j, 1)
670                 qair1(j) = yq(j, 1) + y_d_q(j, 1)
671                 zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
672                      1)))*(ypaprs(j, 1)-ypplay(j, 1))
673                 tairsol(j) = yts(j) + y_d_ts(j)
674                 rugo1(j) = yrugos(j)
675                 IF (nsrf == is_oce) THEN
676                    rugo1(j) = rugos(i, nsrf)
677                 END IF
678                 psfce(j) = ypaprs(j, 1)
679                 patm(j) = ypplay(j, 1)
680    
681         DO j = 1, knon               qairsol(j) = yqsurf(j)
           i = ni(j)  
           DO k = 1, klev  
              d_t(i, k) = d_t(i, k) + y_d_t(j, k)  
              d_q(i, k) = d_q(i, k) + y_d_q(j, k)  
              d_u(i, k) = d_u(i, k) + y_d_u(j, k)  
              d_v(i, k) = d_v(i, k) + y_d_v(j, k)  
              zcoefh(i, k) = zcoefh(i, k) + ycoefh(j, k)  
682            END DO            END DO
        END DO  
   
        !cc diagnostic t, q a 2m et u, v a 10m  
683    
684         DO j = 1, knon            CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &
685            i = ni(j)                 zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &
686            uzon(j) = yu(j, 1) + y_d_u(j, 1)                 yt10m, yq10m, yu10m, yustar)
           vmer(j) = yv(j, 1) + y_d_v(j, 1)  
           tair1(j) = yt(j, 1) + y_d_t(j, 1)  
           qair1(j) = yq(j, 1) + y_d_q(j, 1)  
           zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &  
                1)))*(ypaprs(j, 1)-ypplay(j, 1))  
           tairsol(j) = yts(j) + y_d_ts(j)  
           rugo1(j) = yrugos(j)  
           IF (nsrf == is_oce) THEN  
              rugo1(j) = rugos(i, nsrf)  
           END IF  
           psfce(j) = ypaprs(j, 1)  
           patm(j) = ypplay(j, 1)  
687    
688            qairsol(j) = yqsurf(j)            DO j = 1, knon
689         END DO               i = ni(j)
690                 t2m(i, nsrf) = yt2m(j)
691                 q2m(i, nsrf) = yq2m(j)
692    
693         CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &               ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
694              tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &               u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
695              yu10m, yustar)               v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)
   
        DO j = 1, knon  
           i = ni(j)  
           t2m(i, nsrf) = yt2m(j)  
           q2m(i, nsrf) = yq2m(j)  
   
           ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman  
           u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)  
           v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)  
696    
697         END DO            END DO
698    
699         DO i = 1, knon            CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &
700            y_cd_h(i) = ycoefh(i, 1)                 y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &
701            y_cd_m(i) = ycoefm(i, 1)                 ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
        END DO  
        CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &  
             y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &  
             ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)  
   
        DO j = 1, knon  
           i = ni(j)  
           pblh(i, nsrf) = ypblh(j)  
           plcl(i, nsrf) = ylcl(j)  
           capcl(i, nsrf) = ycapcl(j)  
           oliqcl(i, nsrf) = yoliqcl(j)  
           cteicl(i, nsrf) = ycteicl(j)  
           pblt(i, nsrf) = ypblt(j)  
           therm(i, nsrf) = ytherm(j)  
           trmb1(i, nsrf) = ytrmb1(j)  
           trmb2(i, nsrf) = ytrmb2(j)  
           trmb3(i, nsrf) = ytrmb3(j)  
        END DO  
702    
        DO j = 1, knon  
           DO k = 1, klev + 1  
              i = ni(j)  
              q2(i, k, nsrf) = yq2(j, k)  
           END DO  
        END DO  
        !IM "slab" ocean  
        IF (nsrf == is_oce) THEN  
703            DO j = 1, knon            DO j = 1, knon
              ! on projette sur la grille globale  
704               i = ni(j)               i = ni(j)
705               IF (pctsrf_new(i, is_oce)>epsfra) THEN               pblh(i, nsrf) = ypblh(j)
706                  flux_o(i) = y_flux_o(j)               plcl(i, nsrf) = ylcl(j)
707               ELSE               capcl(i, nsrf) = ycapcl(j)
708                  flux_o(i) = 0.               oliqcl(i, nsrf) = yoliqcl(j)
709               END IF               cteicl(i, nsrf) = ycteicl(j)
710                 pblt(i, nsrf) = ypblt(j)
711                 therm(i, nsrf) = ytherm(j)
712                 trmb1(i, nsrf) = ytrmb1(j)
713                 trmb2(i, nsrf) = ytrmb2(j)
714                 trmb3(i, nsrf) = ytrmb3(j)
715            END DO            END DO
        END IF  
716    
        IF (nsrf == is_sic) THEN  
717            DO j = 1, knon            DO j = 1, knon
718               i = ni(j)               DO k = 1, klev + 1
719               ! On pondère lorsque l'on fait le bilan au sol :                  i = ni(j)
720               IF (pctsrf_new(i, is_sic)>epsfra) THEN                  q2(i, k, nsrf) = yq2(j, k)
721                  flux_g(i) = y_flux_g(j)               END DO
              ELSE  
                 flux_g(i) = 0.  
              END IF  
722            END DO            END DO
723              !IM "slab" ocean
        END IF  
        IF (ocean == 'slab  ') THEN  
724            IF (nsrf == is_oce) THEN            IF (nsrf == is_oce) THEN
725               tslab(1:klon) = ytslab(1:klon)               DO j = 1, knon
726               seaice(1:klon) = y_seaice(1:klon)                  ! on projette sur la grille globale
727                    i = ni(j)
728                    IF (pctsrf_new(i, is_oce)>epsfra) THEN
729                       flux_o(i) = y_flux_o(j)
730                    ELSE
731                       flux_o(i) = 0.
732                    END IF
733                 END DO
734            END IF            END IF
735         END IF  
736              IF (nsrf == is_sic) THEN
737                 DO j = 1, knon
738                    i = ni(j)
739                    ! On pondère lorsque l'on fait le bilan au sol :
740                    IF (pctsrf_new(i, is_sic)>epsfra) THEN
741                       flux_g(i) = y_flux_g(j)
742                    ELSE
743                       flux_g(i) = 0.
744                    END IF
745                 END DO
746    
747              END IF
748              IF (ocean == 'slab  ') THEN
749                 IF (nsrf == is_oce) THEN
750                    tslab(1:klon) = ytslab(1:klon)
751                    seaice(1:klon) = y_seaice(1:klon)
752                 END IF
753              END IF
754           end IF if_knon
755      END DO loop_surface      END DO loop_surface
756    
757      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
Legend:
Removed from v.57  
changed lines
  Added in v.82
  ViewVC Help
Powered by ViewVC 1.1.21