/[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

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

Legend:
Removed from v.47  
changed lines
  Added in v.71

  ViewVC Help
Powered by ViewVC 1.1.21