/[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 47 by guez, Fri Jul 1 15:00:48 2011 UTC trunk/Sources/phylmd/clmain.f revision 154 by guez, Tue Jul 7 17:49:23 2015 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, jour, rmu0, &
8         jour, rmu0, co2_ppm, ok_veget, ocean, npas, nexca, ts,&         co2_ppm, ts, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, &
9         soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil,&         paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat, rain_fall, &
10         qsol, paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat,&         snow_f, solsw, sollw, fder, rlat, rugos, debut, agesno, rugoro, d_t, &
11         rain_f, snow_f, solsw, sollw, sollwdown, fder, rlon, rlat, cufi,&         d_q, d_u, d_v, d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, &
12         cvfi, rugos, debut, lafin, agesno, rugoro, d_t, d_q, d_u, d_v,&         q2, dflux_t, dflux_q, ycoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh, &
13         d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2,&         capcl, oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl, &
14         dflux_t, dflux_q, zcoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh,&         fqcalving, ffonte, run_off_lic_0, flux_o, flux_g, tslab)
15         capcl, oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl,&  
16         fqcalving, ffonte, run_off_lic_0, flux_o, flux_g, tslab, seaice)      ! From phylmd/clmain.F, version 1.6, 2005/11/16 14:47:19
17        ! Author: Z. X. Li (LMD/CNRS), date: 1993/08/18
18      ! From phylmd/clmain.F, version 1.6 2005/11/16 14:47:19      ! Objet : interface de couche limite (diffusion verticale)
19      ! Author: Z.X. Li (LMD/CNRS), date: 1993/08/18  
20      ! Objet : interface de "couche limite" (diffusion verticale)      ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
21        ! de la couche limite pour les traceurs se fait avec "cltrac" et
22      ! Tout ce qui a trait aux traceurs est dans "phytrac" maintenant.      ! ne tient pas compte de la diff\'erentiation des sous-fractions
23      ! Pour l'instant le calcul de la couche limite pour les traceurs      ! de sol.
24      ! se fait avec "cltrac" et ne tient pas compte de la différentiation  
25      ! des sous-fractions de sol.      ! Pour pouvoir extraire les coefficients d'\'echanges et le vent
26        ! dans la premi\`ere couche, trois champs ont \'et\'e cr\'e\'es : "ycoefh",
27      ! Pour pouvoir extraire les coefficients d'échanges et le vent      ! "zu1" et "zv1". Nous avons moyenn\'e les valeurs de ces trois
28      ! dans la première couche, trois champs supplémentaires ont été      ! champs sur les quatre sous-surfaces du mod\`ele.
29      ! créés : "zcoefh", "zu1" et "zv1". Pour l'instant nous avons  
30      ! moyenné les valeurs de ces trois champs sur les 4 sous-surfaces      use clqh_m, only: clqh
31      ! du modèle. Dans l'avenir, si les informations des sous-surfaces      use clvent_m, only: clvent
32      ! doivent être prises en compte, il faudra sortir ces mêmes champs      use coefkz_m, only: coefkz
33      ! en leur ajoutant une dimension, c'est-à-dire "nbsrf" (nombre de      use coefkzmin_m, only: coefkzmin
34      ! sous-surfaces).      USE conf_gcm_m, ONLY: prt_level
35        USE conf_phys_m, ONLY: iflag_pbl
36      ! Arguments:      USE dimens_m, ONLY: iim, jjm
37      ! dtime----input-R- interval du temps (secondes)      USE dimphy, ONLY: klev, klon, zmasq
38      ! itap-----input-I- numero du pas de temps      USE dimsoil, ONLY: nsoilmx
39      ! date0----input-R- jour initial      use hbtm_m, only: hbtm
40      ! t--------input-R- temperature (K)      USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
41      ! q--------input-R- vapeur d'eau (kg/kg)      use stdlevvar_m, only: stdlevvar
42      ! u--------input-R- vitesse u      USE suphec_m, ONLY: rd, rg, rkappa
43      ! v--------input-R- vitesse v      use ustarhb_m, only: ustarhb
44      ! ts-------input-R- temperature du sol (en Kelvin)      use vdif_kcay_m, only: vdif_kcay
45      ! paprs----input-R- pression a intercouche (Pa)      use yamada4_m, only: yamada4
46      ! pplay----input-R- pression au milieu de couche (Pa)  
47      ! radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2      REAL, INTENT(IN):: dtime ! interval du temps (secondes)
48      ! rlat-----input-R- latitude en degree      INTEGER, INTENT(IN):: itap ! numero du pas de temps
49        REAL, INTENT(inout):: pctsrf(klon, nbsrf)
50    
51        ! la nouvelle repartition des surfaces sortie de l'interface
52        REAL, INTENT(out):: pctsrf_new(klon, nbsrf)
53    
54        REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
55        REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg/kg)
56        REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
57        INTEGER, INTENT(IN):: jour ! jour de l'annee en cours
58        REAL, intent(in):: rmu0(klon) ! cosinus de l'angle solaire zenithal    
59        REAL, intent(in):: co2_ppm ! taux CO2 atmosphere
60        REAL, INTENT(IN):: ts(klon, nbsrf) ! temperature du sol (en Kelvin)
61        REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
62        REAL, INTENT(IN):: ksta, ksta_ter
63        LOGICAL, INTENT(IN):: ok_kzmin
64    
65        REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
66        ! soil temperature of surface fraction
67    
68        REAL, INTENT(inout):: qsol(klon)
69        ! column-density of water in soil, in kg m-2
70    
71        REAL, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)
72        REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
73        REAL snow(klon, nbsrf)
74        REAL qsurf(klon, nbsrf)
75        REAL evap(klon, nbsrf)
76        REAL albe(klon, nbsrf)
77        REAL alblw(klon, nbsrf)
78    
79        REAL fluxlat(klon, nbsrf)
80    
81        REAL, intent(in):: rain_fall(klon)
82        ! liquid water mass flux (kg/m2/s), positive down
83    
84        REAL, intent(in):: snow_f(klon)
85        ! solid water mass flux (kg/m2/s), positive down
86    
87        REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)
88        REAL, intent(in):: fder(klon)
89        REAL, INTENT(IN):: rlat(klon) ! latitude en degr\'es
90    
91        REAL rugos(klon, nbsrf)
92      ! rugos----input-R- longeur de rugosite (en m)      ! rugos----input-R- longeur de rugosite (en m)
     ! cufi-----input-R- resolution des mailles en x (m)  
     ! cvfi-----input-R- resolution des mailles en y (m)  
93    
94        LOGICAL, INTENT(IN):: debut
95        real agesno(klon, nbsrf)
96        REAL, INTENT(IN):: rugoro(klon)
97    
98        REAL d_t(klon, klev), d_q(klon, klev)
99      ! d_t------output-R- le changement pour "t"      ! d_t------output-R- le changement pour "t"
100      ! d_q------output-R- le changement pour "q"      ! d_q------output-R- le changement pour "q"
101      ! d_u------output-R- le changement pour "u"  
102      ! d_v------output-R- le changement pour "v"      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
103      ! d_ts-----output-R- le changement pour "ts"      ! changement pour "u" et "v"
104    
105        REAL, intent(out):: d_ts(klon, nbsrf) ! le changement pour "ts"
106    
107        REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)
108      ! 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)
109      !                    (orientation positive vers le bas)      !                    (orientation positive vers le bas)
110      ! 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)
111    
112        REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)
113      ! 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
114      ! 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
115    
116        REAL, INTENT(out):: cdragh(klon), cdragm(klon)
117        real q2(klon, klev+1, nbsrf)
118    
119        REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)
120      ! dflux_t derive du flux sensible      ! dflux_t derive du flux sensible
121      ! dflux_q derive du flux latent      ! dflux_q derive du flux latent
122      !IM "slab" ocean      !IM "slab" ocean
     ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')  
     ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')  
123    
124      ! tslab-in/output-R temperature du slab ocean (en Kelvin)      REAL, intent(out):: ycoefh(klon, klev)
125      ! uniqmnt pour slab      REAL, intent(out):: zu1(klon)
126        REAL zv1(klon)
127        REAL t2m(klon, nbsrf), q2m(klon, nbsrf)
128        REAL u10m(klon, nbsrf), v10m(klon, nbsrf)
129    
130      ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')      !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds
131      !cc      ! physiq ce qui permet de sortir les grdeurs par sous surface)
132      ! ffonte----Flux thermique utilise pour fondre la neige      REAL pblh(klon, nbsrf)
133      ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la      ! pblh------- HCL
134      !           hauteur de neige, en kg/m2/s      REAL capcl(klon, nbsrf)
135      ! on rajoute en output yu1 et yv1 qui sont les vents dans      REAL oliqcl(klon, nbsrf)
136      ! la premiere couche      REAL cteicl(klon, nbsrf)
137      ! ces 4 variables sont maintenant traites dans phytrac      REAL pblt(klon, nbsrf)
138      ! itr--------input-I- nombre de traceurs      ! pblT------- T au nveau HCL
139      ! tr---------input-R- q. de traceurs      REAL therm(klon, nbsrf)
140      ! 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  
141      ! trmb1-------deep_cape      ! trmb1-------deep_cape
142        REAL trmb2(klon, nbsrf)
143      ! trmb2--------inhibition      ! trmb2--------inhibition
144        REAL trmb3(klon, nbsrf)
145      ! trmb3-------Point Omega      ! trmb3-------Point Omega
146      ! 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)  
147      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
148      REAL run_off_lic_0(klon), y_run_off_lic_0(klon)      ! ffonte----Flux thermique utilise pour fondre la neige
149        ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
150      REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)      !           hauteur de neige, en kg/m2/s
151      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)  
   
     REAL fluxlat(klon, nbsrf)  
152    
153      REAL rain_f(klon), snow_f(klon)      REAL flux_o(klon), flux_g(klon)
154      REAL fder(klon)      !IM "slab" ocean
155        ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')
156        ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')
157    
158      REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)      REAL tslab(klon)
159      REAL rugos(klon, nbsrf)      ! tslab-in/output-R temperature du slab ocean (en Kelvin)
160      ! la nouvelle repartition des surfaces sortie de l'interface      ! uniqmnt pour slab
     REAL pctsrf_new(klon, nbsrf)  
161    
162      REAL zcoefh(klon, klev)      ! Local:
     REAL zu1(klon)  
     REAL zv1(klon)  
163    
164      !$$$ PB ajout pour soil      REAL y_flux_o(klon), y_flux_g(klon)
165      LOGICAL, INTENT (IN) :: soil_model      real ytslab(klon)
166      !IM ajout seuils cdrm, cdrh      REAL y_fqcalving(klon), y_ffonte(klon)
167      REAL cdmmax, cdhmax      real y_run_off_lic_0(klon)
168    
169      REAL ksta, ksta_ter      REAL rugmer(klon)
     LOGICAL ok_kzmin  
170    
     REAL ftsoil(klon, nsoilmx, nbsrf)  
171      REAL ytsoil(klon, nsoilmx)      REAL ytsoil(klon, nsoilmx)
     REAL qsol(klon)  
   
     EXTERNAL clqh, clvent, calbeta, cltrac  
172    
173      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
174      REAL yalb(klon)      REAL yalb(klon)
175      REAL yalblw(klon)      REAL yalblw(klon)
176      REAL yu1(klon), yv1(klon)      REAL yu1(klon), yv1(klon)
177      REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)      ! on rajoute en output yu1 et yv1 qui sont les vents dans
178      REAL yrain_f(klon), ysnow_f(klon)      ! la premiere couche
179      REAL ysollw(klon), ysolsw(klon), ysollwdown(klon)      REAL ysnow(klon), yqsurf(klon), yagesno(klon)
180      REAL yfder(klon), ytaux(klon), ytauy(klon)  
181        real yqsol(klon)
182        ! column-density of water in soil, in kg m-2
183    
184        REAL yrain_f(klon)
185        ! liquid water mass flux (kg/m2/s), positive down
186    
187        REAL ysnow_f(klon)
188        ! solid water mass flux (kg/m2/s), positive down
189    
190        REAL ysollw(klon), ysolsw(klon)
191        REAL yfder(klon)
192      REAL yrugm(klon), yrads(klon), yrugoro(klon)      REAL yrugm(klon), yrads(klon), yrugoro(klon)
193    
194      REAL yfluxlat(klon)      REAL yfluxlat(klon)
# Line 202  contains Line 199  contains
199      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)
200      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)
201      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
202      REAL ycoefh(klon, klev), ycoefm(klon, klev)      REAL coefh(klon, klev), coefm(klon, klev)
203      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
204      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
205      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)
206    
     LOGICAL ok_nonloc  
     PARAMETER (ok_nonloc=.FALSE.)  
207      REAL ycoefm0(klon, klev), ycoefh0(klon, klev)      REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
208    
209      REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)      REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
210      REAL ykmm(klon, klev+1), ykmn(klon, klev+1)      REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
211      REAL ykmq(klon, klev+1)      REAL ykmq(klon, klev+1)
212      REAL yq2(klon, klev+1), q2(klon, klev+1, nbsrf)      REAL yq2(klon, klev+1)
213      REAL q2diag(klon, klev+1)      REAL q2diag(klon, klev+1)
214    
215      REAL u1lay(klon), v1lay(klon)      REAL u1lay(klon), v1lay(klon)
# Line 224  contains Line 219  contains
219      INTEGER ni(klon), knon, j      INTEGER ni(klon), knon, j
220    
221      REAL pctsrf_pot(klon, nbsrf)      REAL pctsrf_pot(klon, nbsrf)
222      ! "pourcentage potentiel" pour tenir compte des éventuelles      ! "pourcentage potentiel" pour tenir compte des \'eventuelles
223      ! apparitions ou disparitions de la glace de mer      ! apparitions ou disparitions de la glace de mer
224    
225      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
226    
     ! maf pour sorties IOISPL en cas de debugagage  
   
     CHARACTER (80) cldebug  
     SAVE cldebug  
     CHARACTER (8) cl_surf(nbsrf)  
     SAVE cl_surf  
     INTEGER nhoridbg, nidbg  
     SAVE nhoridbg, nidbg  
     INTEGER ndexbg(iim*(jjm+1))  
     REAL zx_lon(iim, jjm+1), zx_lat(iim, jjm+1), zjulian  
     REAL tabindx(klon)  
     REAL debugtab(iim, jjm+1)  
     LOGICAL first_appel  
     SAVE first_appel  
     DATA first_appel/ .TRUE./  
     LOGICAL :: debugindex = .FALSE.  
     INTEGER idayref  
     REAL t2m(klon, nbsrf), q2m(klon, nbsrf)  
     REAL u10m(klon, nbsrf), v10m(klon, nbsrf)  
   
227      REAL yt2m(klon), yq2m(klon), yu10m(klon)      REAL yt2m(klon), yq2m(klon), yu10m(klon)
228      REAL yustar(klon)      REAL yustar(klon)
229      ! -- LOOP      ! -- LOOP
# Line 258  contains Line 233  contains
233      ! -- LOOP      ! -- LOOP
234    
235      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)  
236      REAL ypblh(klon)      REAL ypblh(klon)
237      REAL ylcl(klon)      REAL ylcl(klon)
238      REAL ycapcl(klon)      REAL ycapcl(klon)
# Line 280  contains Line 243  contains
243      REAL ytrmb1(klon)      REAL ytrmb1(klon)
244      REAL ytrmb2(klon)      REAL ytrmb2(klon)
245      REAL ytrmb3(klon)      REAL ytrmb3(klon)
     REAL y_cd_h(klon), y_cd_m(klon)  
246      REAL uzon(klon), vmer(klon)      REAL uzon(klon), vmer(klon)
247      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
248      REAL psfce(klon), patm(klon)      REAL psfce(klon), patm(klon)
# Line 292  contains Line 254  contains
254      LOGICAL zxli      LOGICAL zxli
255      PARAMETER (zxli=.FALSE.)      PARAMETER (zxli=.FALSE.)
256    
     REAL zt, zqs, zdelta, zcor  
     REAL t_coup  
     PARAMETER (t_coup=273.15)  
   
     CHARACTER (len=20) :: modname = 'clmain'  
   
257      !------------------------------------------------------------      !------------------------------------------------------------
258    
259      ytherm = 0.      ytherm = 0.
260    
     IF (debugindex .AND. first_appel) THEN  
        first_appel = .FALSE.  
   
        ! initialisation sorties netcdf  
   
        idayref = day_ini  
        CALL ymds2ju(annee_ref, 1, idayref, 0., zjulian)  
        CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlon, zx_lon)  
        DO i = 1, iim  
           zx_lon(i, 1) = rlon(i+1)  
           zx_lon(i, jjm+1) = rlon(i+1)  
        END DO  
        CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlat, zx_lat)  
        cldebug = 'sous_index'  
        CALL histbeg_totreg(cldebug, zx_lon(:, 1), zx_lat(1, :), 1, &  
             iim, 1, jjm+1, itau_phy, zjulian, dtime, nhoridbg, nidbg)  
        ! no vertical axis  
        cl_surf(1) = 'ter'  
        cl_surf(2) = 'lic'  
        cl_surf(3) = 'oce'  
        cl_surf(4) = 'sic'  
        DO nsrf = 1, nbsrf  
           CALL histdef(nidbg, cl_surf(nsrf), cl_surf(nsrf), '-', iim, jjm+1, &  
                nhoridbg, 1, 1, 1, -99, 'inst', dtime, dtime)  
        END DO  
        CALL histend(nidbg)  
        CALL histsync(nidbg)  
     END IF  
   
261      DO k = 1, klev ! epaisseur de couche      DO k = 1, klev ! epaisseur de couche
262         DO i = 1, klon         DO i = 1, klon
263            delp(i, k) = paprs(i, k) - paprs(i, k+1)            delp(i, k) = paprs(i, k) - paprs(i, k+1)
# Line 356  contains Line 283  contains
283      ysnow = 0.      ysnow = 0.
284      yqsurf = 0.      yqsurf = 0.
285      yalb = 0.      yalb = 0.
     yalblw = 0.  
286      yrain_f = 0.      yrain_f = 0.
287      ysnow_f = 0.      ysnow_f = 0.
288      yfder = 0.      yfder = 0.
     ytaux = 0.  
     ytauy = 0.  
289      ysolsw = 0.      ysolsw = 0.
290      ysollw = 0.      ysollw = 0.
     ysollwdown = 0.  
291      yrugos = 0.      yrugos = 0.
292      yu1 = 0.      yu1 = 0.
293      yv1 = 0.      yv1 = 0.
# Line 379  contains Line 302  contains
302      pctsrf_new = 0.      pctsrf_new = 0.
303      y_flux_u = 0.      y_flux_u = 0.
304      y_flux_v = 0.      y_flux_v = 0.
     !$$ PB  
305      y_dflux_t = 0.      y_dflux_t = 0.
306      y_dflux_q = 0.      y_dflux_q = 0.
307      ytsoil = 999999.      ytsoil = 999999.
308      yrugoro = 0.      yrugoro = 0.
     ! -- LOOP  
309      yu10mx = 0.      yu10mx = 0.
310      yu10my = 0.      yu10my = 0.
311      ywindsp = 0.      ywindsp = 0.
     ! -- LOOP  
312      d_ts = 0.      d_ts = 0.
     !§§§ PB  
313      yfluxlat = 0.      yfluxlat = 0.
314      flux_t = 0.      flux_t = 0.
315      flux_q = 0.      flux_q = 0.
# Line 400  contains Line 319  contains
319      d_q = 0.      d_q = 0.
320      d_u = 0.      d_u = 0.
321      d_v = 0.      d_v = 0.
322      zcoefh = 0.      ycoefh = 0.
   
     ! Boucler sur toutes les sous-fractions du sol:  
323    
324      ! Initialisation des "pourcentages potentiels". On considère ici qu'on      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
325      ! peut avoir potentiellement de la glace sur tout le domaine océanique      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
326      ! (à affiner)      ! (\`a affiner)
327    
328      pctsrf_pot = pctsrf      pctsrf_pot = pctsrf
329      pctsrf_pot(:, is_oce) = 1. - zmasq      pctsrf_pot(:, is_oce) = 1. - zmasq
330      pctsrf_pot(:, is_sic) = 1. - zmasq      pctsrf_pot(:, is_sic) = 1. - zmasq
331    
332      DO nsrf = 1, nbsrf      ! Boucler sur toutes les sous-fractions du sol:
333         ! chercher les indices:  
334        loop_surface: DO nsrf = 1, nbsrf
335           ! Chercher les indices :
336         ni = 0         ni = 0
337         knon = 0         knon = 0
338         DO i = 1, klon         DO i = 1, klon
339            ! Pour déterminer le domaine à traiter, on utilise les surfaces            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
340            ! "potentielles"            ! "potentielles"
341            IF (pctsrf_pot(i, nsrf) > epsfra) THEN            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
342               knon = knon + 1               knon = knon + 1
# Line 425  contains Line 344  contains
344            END IF            END IF
345         END DO         END DO
346    
347         ! variables pour avoir une sortie IOIPSL des INDEX         if_knon: IF (knon /= 0) then
        IF (debugindex) THEN  
           tabindx = 0.  
           DO i = 1, knon  
              tabindx(i) = real(i)  
           END DO  
           debugtab = 0.  
           ndexbg = 0  
           CALL gath2cpl(tabindx, debugtab, klon, knon, iim, jjm, ni)  
           CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)  
        END IF  
   
        IF (knon == 0) CYCLE  
   
        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  
348            DO j = 1, knon            DO j = 1, knon
349               i = ni(j)               i = ni(j)
350               ypaprs(j, k) = paprs(i, k)               ypct(j) = pctsrf(i, nsrf)
351               ypplay(j, k) = pplay(i, k)               yts(j) = ts(i, nsrf)
352               ydelp(j, k) = delp(i, k)               ytslab(i) = tslab(i)
353               yu(j, k) = u(i, k)               ysnow(j) = snow(i, nsrf)
354               yv(j, k) = v(i, k)               yqsurf(j) = qsurf(i, nsrf)
355               yt(j, k) = t(i, k)               yalb(j) = albe(i, nsrf)
356               yq(j, k) = q(i, k)               yrain_f(j) = rain_fall(i)
357            END DO               ysnow_f(j) = snow_f(i)
358         END DO               yagesno(j) = agesno(i, nsrf)
359                 yfder(j) = fder(i)
360                 ysolsw(j) = solsw(i, nsrf)
361                 ysollw(j) = sollw(i, nsrf)
362                 yrugos(j) = rugos(i, nsrf)
363                 yrugoro(j) = rugoro(i)
364                 yu1(j) = u1lay(i)
365                 yv1(j) = v1lay(i)
366                 yrads(j) = ysolsw(j) + ysollw(j)
367                 ypaprs(j, klev+1) = paprs(i, klev+1)
368                 y_run_off_lic_0(j) = run_off_lic_0(i)
369                 yu10mx(j) = u10m(i, nsrf)
370                 yu10my(j) = v10m(i, nsrf)
371                 ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))
372              END DO
373    
374              ! For continent, copy soil water content
375              IF (nsrf == is_ter) THEN
376                 yqsol(:knon) = qsol(ni(:knon))
377              ELSE
378                 yqsol = 0.
379              END IF
380    
381         ! calculer Cdrag et les coefficients d'echange            DO k = 1, nsoilmx
382         CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,&               DO j = 1, knon
383              yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)                  i = ni(j)
384         IF (iflag_pbl == 1) THEN                  ytsoil(j, k) = ftsoil(i, k, nsrf)
           CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)  
           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))  
385               END DO               END DO
386            END DO            END DO
        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)  
387    
388            DO k = 1, klev            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))  
              END DO  
           END DO  
        END IF  
   
        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  
           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  
389               DO j = 1, knon               DO j = 1, knon
390                  i = ni(j)                  i = ni(j)
391                  yq2(j, k) = q2(i, k, nsrf)                  ypaprs(j, k) = paprs(i, k)
392                    ypplay(j, k) = pplay(i, k)
393                    ydelp(j, k) = delp(i, k)
394                    yu(j, k) = u(i, k)
395                    yv(j, k) = v(i, k)
396                    yt(j, k) = t(i, k)
397                    yq(j, k) = q(i, k)
398               END DO               END DO
399            END DO            END DO
400    
401            y_cd_m(1:knon) = ycoefm(1:knon, 1)            ! calculer Cdrag et les coefficients d'echange
402            y_cd_h(1:knon) = ycoefh(1:knon, 1)            CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &
403            CALL ustarhb(knon, yu, yv, y_cd_m, yustar)                 yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))
404              IF (iflag_pbl == 1) THEN
405                 CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
406                 coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
407                 coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
408              END IF
409    
410            IF (prt_level>9) THEN            ! on met un seuil pour coefm et coefh
411               PRINT *, 'USTAR = ', yustar            IF (nsrf == is_oce) THEN
412                 coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
413                 coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
414            END IF            END IF
415    
416            ! iflag_pbl peut être utilisé comme longueur de mélange            IF (ok_kzmin) THEN
417                 ! Calcul d'une diffusion minimale pour les conditions tres stables
418                 CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
419                      coefm(:knon, 1), ycoefm0, ycoefh0)
420                 coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
421                 coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
422              END IF
423    
424            IF (iflag_pbl >= 11) THEN            IF (iflag_pbl >= 3) THEN
425               CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &               ! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et
426                    yu, yv, yteta, y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, &               ! Fr\'ed\'eric Hourdin
427                    iflag_pbl)               yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
428            ELSE                    + ypplay(:knon, 1))) &
429               CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &                    * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
430                    y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)               DO k = 2, klev
431                    yzlay(1:knon, k) = yzlay(1:knon, k-1) &
432                         + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
433                         / ypaprs(1:knon, k) &
434                         * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
435                 END DO
436                 DO k = 1, klev
437                    yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
438                         / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
439                 END DO
440                 yzlev(1:knon, 1) = 0.
441                 yzlev(:knon, klev+1) = 2. * yzlay(:knon, klev) &
442                      - yzlay(:knon, klev - 1)
443                 DO k = 2, klev
444                    yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
445                 END DO
446                 DO k = 1, klev + 1
447                    DO j = 1, knon
448                       i = ni(j)
449                       yq2(j, k) = q2(i, k, nsrf)
450                    END DO
451                 END DO
452    
453                 CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
454                 IF (prt_level > 9) PRINT *, 'USTAR = ', yustar
455    
456                 ! iflag_pbl peut \^etre utilis\'e comme longueur de m\'elange
457    
458                 IF (iflag_pbl >= 11) THEN
459                    CALL vdif_kcay(knon, dtime, rg, ypaprs, yzlev, yzlay, yu, yv, &
460                         yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, yustar, &
461                         iflag_pbl)
462                 ELSE
463                    CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
464                         coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
465                 END IF
466    
467                 coefm(:knon, 2:) = ykmm(:knon, 2:klev)
468                 coefh(:knon, 2:) = ykmn(:knon, 2:klev)
469            END IF            END IF
470    
471            ycoefm(1:knon, 1) = y_cd_m(1:knon)            ! calculer la diffusion des vitesses "u" et "v"
472            ycoefh(1:knon, 1) = y_cd_h(1:knon)            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &
473            ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)                 ypplay, ydelp, y_d_u, y_flux_u)
474            ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
475         END IF                 ypplay, ydelp, y_d_v, y_flux_v)
476    
477         ! calculer la diffusion des vitesses "u" et "v"            ! calculer la diffusion de "q" et de "h"
478         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &            CALL clqh(dtime, itap, jour, debut, rlat, knon, nsrf, ni(:knon), &
479              ydelp, y_d_u, y_flux_u)                 pctsrf, ytsoil, yqsol, rmu0, co2_ppm, yrugos, yrugoro, yu1, &
480         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &                 yv1, coefh(:knon, :), yt, yq, yts, ypaprs, ypplay, ydelp, &
481              ydelp, y_d_v, y_flux_v)                 yrads, yalb, yalblw(:knon), ysnow, yqsurf, yrain_f, ysnow_f, &
482                   yfder, ysolsw, yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, &
483         ! pour le couplage                 y_d_ts(:knon), yz0_new, y_flux_t, y_flux_q, y_dflux_t, &
484         ytaux = y_flux_u(:, 1)                 y_dflux_q, y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, &
485         ytauy = y_flux_v(:, 1)                 y_flux_g)
   
        ! 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  
           DO j = 1, knon  
              yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &  
                   0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2))  
              yrugm(j) = max(1.5E-05, yrugm(j))  
           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  
486    
487         DO k = 1, klev            ! calculer la longueur de rugosite sur ocean
488              yrugm = 0.
489              IF (nsrf == is_oce) THEN
490                 DO j = 1, knon
491                    yrugm(j) = 0.018*coefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
492                         0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))
493                    yrugm(j) = max(1.5E-05, yrugm(j))
494                 END DO
495              END IF
496            DO j = 1, knon            DO j = 1, knon
497               i = ni(j)               y_dflux_t(j) = y_dflux_t(j)*ypct(j)
498               ycoefh(j, k) = ycoefh(j, k)*ypct(j)               y_dflux_q(j) = y_dflux_q(j)*ypct(j)
499               ycoefm(j, k) = ycoefm(j, k)*ypct(j)               yu1(j) = yu1(j)*ypct(j)
500               y_d_t(j, k) = y_d_t(j, k)*ypct(j)               yv1(j) = yv1(j)*ypct(j)
              y_d_q(j, k) = y_d_q(j, k)*ypct(j)  
              flux_t(i, k, nsrf) = y_flux_t(j, k)  
              flux_q(i, k, nsrf) = y_flux_q(j, k)  
              flux_u(i, k, nsrf) = y_flux_u(j, k)  
              flux_v(i, k, nsrf) = y_flux_v(j, k)  
              y_d_u(j, k) = y_d_u(j, k)*ypct(j)  
              y_d_v(j, k) = y_d_v(j, k)*ypct(j)  
501            END DO            END DO
        END DO  
502    
503         evap(:, nsrf) = -flux_q(:, 1, nsrf)            DO k = 1, klev
504                 DO j = 1, knon
505                    i = ni(j)
506                    coefh(j, k) = coefh(j, k)*ypct(j)
507                    coefm(j, k) = coefm(j, k)*ypct(j)
508                    y_d_t(j, k) = y_d_t(j, k)*ypct(j)
509                    y_d_q(j, k) = y_d_q(j, k)*ypct(j)
510                    flux_t(i, k, nsrf) = y_flux_t(j, k)
511                    flux_q(i, k, nsrf) = y_flux_q(j, k)
512                    flux_u(i, k, nsrf) = y_flux_u(j, k)
513                    flux_v(i, k, nsrf) = y_flux_v(j, k)
514                    y_d_u(j, k) = y_d_u(j, k)*ypct(j)
515                    y_d_v(j, k) = y_d_v(j, k)*ypct(j)
516                 END DO
517              END DO
518    
519         albe(:, nsrf) = 0.            evap(:, nsrf) = -flux_q(:, 1, nsrf)
520         alblw(:, nsrf) = 0.  
521         snow(:, nsrf) = 0.            albe(:, nsrf) = 0.
522         qsurf(:, nsrf) = 0.            alblw(:, nsrf) = 0.
523         rugos(:, nsrf) = 0.            snow(:, nsrf) = 0.
524         fluxlat(:, nsrf) = 0.            qsurf(:, nsrf) = 0.
525         DO j = 1, knon            rugos(:, nsrf) = 0.
526            i = ni(j)            fluxlat(:, nsrf) = 0.
           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  
527            DO j = 1, knon            DO j = 1, knon
528               i = ni(j)               i = ni(j)
529               qsol(i) = yqsol(j)               d_ts(i, nsrf) = y_d_ts(j)
530                 albe(i, nsrf) = yalb(j)
531                 alblw(i, nsrf) = yalblw(j)
532                 snow(i, nsrf) = ysnow(j)
533                 qsurf(i, nsrf) = yqsurf(j)
534                 rugos(i, nsrf) = yz0_new(j)
535                 fluxlat(i, nsrf) = yfluxlat(j)
536                 IF (nsrf == is_oce) THEN
537                    rugmer(i) = yrugm(j)
538                    rugos(i, nsrf) = yrugm(j)
539                 END IF
540                 agesno(i, nsrf) = yagesno(j)
541                 fqcalving(i, nsrf) = y_fqcalving(j)
542                 ffonte(i, nsrf) = y_ffonte(j)
543                 cdragh(i) = cdragh(i) + coefh(j, 1)
544                 cdragm(i) = cdragm(i) + coefm(j, 1)
545                 dflux_t(i) = dflux_t(i) + y_dflux_t(j)
546                 dflux_q(i) = dflux_q(i) + y_dflux_q(j)
547                 zu1(i) = zu1(i) + yu1(j)
548                 zv1(i) = zv1(i) + yv1(j)
549              END DO
550              IF (nsrf == is_ter) THEN
551                 qsol(ni(:knon)) = yqsol(:knon)
552              else IF (nsrf == is_lic) THEN
553                 DO j = 1, knon
554                    i = ni(j)
555                    run_off_lic_0(i) = y_run_off_lic_0(j)
556                 END DO
557              END IF
558    
559              ftsoil(:, :, nsrf) = 0.
560              DO k = 1, nsoilmx
561                 DO j = 1, knon
562                    i = ni(j)
563                    ftsoil(i, k, nsrf) = ytsoil(j, k)
564                 END DO
565            END DO            END DO
566         END IF  
        IF (nsrf == is_lic) THEN  
567            DO j = 1, knon            DO j = 1, knon
568               i = ni(j)               i = ni(j)
569               run_off_lic_0(i) = y_run_off_lic_0(j)               DO k = 1, klev
570                    d_t(i, k) = d_t(i, k) + y_d_t(j, k)
571                    d_q(i, k) = d_q(i, k) + y_d_q(j, k)
572                    d_u(i, k) = d_u(i, k) + y_d_u(j, k)
573                    d_v(i, k) = d_v(i, k) + y_d_v(j, k)
574                    ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
575                 END DO
576            END DO            END DO
577         END IF  
578         !$$$ PB ajout pour soil            ! diagnostic t, q a 2m et u, v a 10m
579         ftsoil(:, :, nsrf) = 0.  
        DO k = 1, nsoilmx  
580            DO j = 1, knon            DO j = 1, knon
581               i = ni(j)               i = ni(j)
582               ftsoil(i, k, nsrf) = ytsoil(j, k)               uzon(j) = yu(j, 1) + y_d_u(j, 1)
583            END DO               vmer(j) = yv(j, 1) + y_d_v(j, 1)
584         END DO               tair1(j) = yt(j, 1) + y_d_t(j, 1)
585                 qair1(j) = yq(j, 1) + y_d_q(j, 1)
586                 zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
587                      1)))*(ypaprs(j, 1)-ypplay(j, 1))
588                 tairsol(j) = yts(j) + y_d_ts(j)
589                 rugo1(j) = yrugos(j)
590                 IF (nsrf == is_oce) THEN
591                    rugo1(j) = rugos(i, nsrf)
592                 END IF
593                 psfce(j) = ypaprs(j, 1)
594                 patm(j) = ypplay(j, 1)
595    
596         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)  
597            END DO            END DO
        END DO  
598    
599         !cc diagnostic t, q a 2m et u, v a 10m            CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &
600                   zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &
601                   yt10m, yq10m, yu10m, yustar)
602    
603         DO j = 1, knon            DO j = 1, knon
604            i = ni(j)               i = ni(j)
605            uzon(j) = yu(j, 1) + y_d_u(j, 1)               t2m(i, nsrf) = yt2m(j)
606            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  
607    
608         CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &               ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
609              tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &               u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
610              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)  
611    
612         END DO            END DO
613    
614         DO i = 1, knon            CALL hbtm(knon, ypaprs, ypplay, yt2m, yq2m, yustar, &
615            y_cd_h(i) = ycoefh(i, 1)                 y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &
616            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  
617    
        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  
618            DO j = 1, knon            DO j = 1, knon
              ! on projette sur la grille globale  
619               i = ni(j)               i = ni(j)
620               IF (pctsrf_new(i, is_oce)>epsfra) THEN               pblh(i, nsrf) = ypblh(j)
621                  flux_o(i) = y_flux_o(j)               plcl(i, nsrf) = ylcl(j)
622               ELSE               capcl(i, nsrf) = ycapcl(j)
623                  flux_o(i) = 0.               oliqcl(i, nsrf) = yoliqcl(j)
624               END IF               cteicl(i, nsrf) = ycteicl(j)
625                 pblt(i, nsrf) = ypblt(j)
626                 therm(i, nsrf) = ytherm(j)
627                 trmb1(i, nsrf) = ytrmb1(j)
628                 trmb2(i, nsrf) = ytrmb2(j)
629                 trmb3(i, nsrf) = ytrmb3(j)
630            END DO            END DO
        END IF  
631    
        IF (nsrf == is_sic) THEN  
632            DO j = 1, knon            DO j = 1, knon
633               i = ni(j)               DO k = 1, klev + 1
634               ! On pondère lorsque l'on fait le bilan au sol :                  i = ni(j)
635               IF (pctsrf_new(i, is_sic)>epsfra) THEN                  q2(i, k, nsrf) = yq2(j, k)
636                  flux_g(i) = y_flux_g(j)               END DO
              ELSE  
                 flux_g(i) = 0.  
              END IF  
637            END DO            END DO
638              !IM "slab" ocean
        END IF  
        IF (ocean == 'slab  ') THEN  
639            IF (nsrf == is_oce) THEN            IF (nsrf == is_oce) THEN
640               tslab(1:klon) = ytslab(1:klon)               DO j = 1, knon
641               seaice(1:klon) = y_seaice(1:klon)                  ! on projette sur la grille globale
642                    i = ni(j)
643                    IF (pctsrf_new(i, is_oce)>epsfra) THEN
644                       flux_o(i) = y_flux_o(j)
645                    ELSE
646                       flux_o(i) = 0.
647                    END IF
648                 END DO
649            END IF            END IF
650         END IF  
651      END DO            IF (nsrf == is_sic) THEN
652                 DO j = 1, knon
653                    i = ni(j)
654                    ! On pond\`ere lorsque l'on fait le bilan au sol :
655                    IF (pctsrf_new(i, is_sic)>epsfra) THEN
656                       flux_g(i) = y_flux_g(j)
657                    ELSE
658                       flux_g(i) = 0.
659                    END IF
660                 END DO
661    
662              END IF
663           end IF if_knon
664        END DO loop_surface
665    
666      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
667    

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

  ViewVC Help
Powered by ViewVC 1.1.21