/[lmdze]/trunk/phylmd/pbl_surface.f
ViewVC logotype

Diff of /trunk/phylmd/pbl_surface.f

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

trunk/libf/phylmd/clmain.f90 revision 40 by guez, Tue Feb 22 13:49:36 2011 UTC trunk/Sources/phylmd/clmain.f revision 155 by guez, Wed Jul 8 17:03:45 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, falbe, fluxlat, rain_fall, snow_f, &
10         qsol, paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat,&         solsw, sollw, fder, rlat, rugos, debut, agesno, rugoro, d_t, d_q, d_u, &
11         rain_f, snow_f, solsw, sollw, sollwdown, fder, rlon, rlat, cufi,&         d_v, d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, &
12         cvfi, rugos, debut, lafin, agesno, rugoro, d_t, d_q, d_u, d_v,&         dflux_t, dflux_q, ycoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh, capcl, &
13         d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2,&         oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl, fqcalving, &
14         dflux_t, dflux_q, zcoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh,&         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, intent(inout):: falbe(klon, nbsrf)
77    
78        REAL fluxlat(klon, nbsrf)
79    
80        REAL, intent(in):: rain_fall(klon)
81        ! liquid water mass flux (kg/m2/s), positive down
82    
83        REAL, intent(in):: snow_f(klon)
84        ! solid water mass flux (kg/m2/s), positive down
85    
86        REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)
87        REAL, intent(in):: fder(klon)
88        REAL, INTENT(IN):: rlat(klon) ! latitude en degr\'es
89    
90        REAL rugos(klon, nbsrf)
91      ! 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)  
92    
93        LOGICAL, INTENT(IN):: debut
94        real agesno(klon, nbsrf)
95        REAL, INTENT(IN):: rugoro(klon)
96    
97        REAL d_t(klon, klev), d_q(klon, klev)
98      ! d_t------output-R- le changement pour "t"      ! d_t------output-R- le changement pour "t"
99      ! d_q------output-R- le changement pour "q"      ! d_q------output-R- le changement pour "q"
100      ! d_u------output-R- le changement pour "u"  
101      ! d_v------output-R- le changement pour "v"      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
102      ! d_ts-----output-R- le changement pour "ts"      ! changement pour "u" et "v"
103    
104        REAL, intent(out):: d_ts(klon, nbsrf) ! le changement pour "ts"
105    
106        REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)
107      ! 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)
108      !                    (orientation positive vers le bas)      !                    (orientation positive vers le bas)
109      ! 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)
110    
111        REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)
112      ! 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
113      ! 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
114    
115        REAL, INTENT(out):: cdragh(klon), cdragm(klon)
116        real q2(klon, klev+1, nbsrf)
117    
118        REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)
119      ! dflux_t derive du flux sensible      ! dflux_t derive du flux sensible
120      ! dflux_q derive du flux latent      ! dflux_q derive du flux latent
121      !IM "slab" ocean      !IM "slab" ocean
     ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')  
     ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')  
122    
123      ! tslab-in/output-R temperature du slab ocean (en Kelvin)      REAL, intent(out):: ycoefh(klon, klev)
124      ! uniqmnt pour slab      REAL, intent(out):: zu1(klon)
125        REAL zv1(klon)
126        REAL t2m(klon, nbsrf), q2m(klon, nbsrf)
127        REAL u10m(klon, nbsrf), v10m(klon, nbsrf)
128    
129      ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')      !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds
130      !cc      ! physiq ce qui permet de sortir les grdeurs par sous surface)
131      ! ffonte----Flux thermique utilise pour fondre la neige      REAL pblh(klon, nbsrf)
132      ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la      ! pblh------- HCL
133      !           hauteur de neige, en kg/m2/s      REAL capcl(klon, nbsrf)
134      ! on rajoute en output yu1 et yv1 qui sont les vents dans      REAL oliqcl(klon, nbsrf)
135      ! la premiere couche      REAL cteicl(klon, nbsrf)
136      ! ces 4 variables sont maintenant traites dans phytrac      REAL pblt(klon, nbsrf)
137      ! itr--------input-I- nombre de traceurs      ! pblT------- T au nveau HCL
138      ! tr---------input-R- q. de traceurs      REAL therm(klon, nbsrf)
139      ! 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  
140      ! trmb1-------deep_cape      ! trmb1-------deep_cape
141        REAL trmb2(klon, nbsrf)
142      ! trmb2--------inhibition      ! trmb2--------inhibition
143        REAL trmb3(klon, nbsrf)
144      ! trmb3-------Point Omega      ! trmb3-------Point Omega
145      ! 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 histcom, ONLY : histbeg_totreg, histdef, histend, histsync  
     use histwrite_m, only: histwrite  
     use calendar, ONLY : ymds2ju  
     USE dimens_m, ONLY : iim, jjm  
     USE indicesol, ONLY : epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf  
     USE dimphy, ONLY : klev, klon, zmasq  
     USE dimsoil, ONLY : nsoilmx  
     USE temps, ONLY : annee_ref, itau_phy  
     USE dynetat0_m, ONLY : day_ini  
     USE iniprint, ONLY : prt_level  
     USE suphec_m, ONLY : rd, rg, rkappa  
     USE conf_phys_m, ONLY : iflag_pbl  
     USE gath_cpl, ONLY : gath2cpl  
     use hbtm_m, only: hbtm  
   
     REAL, INTENT (IN) :: dtime  
     REAL date0  
     INTEGER, INTENT (IN) :: itap  
     REAL t(klon, klev), q(klon, klev)  
     REAL 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)  
146      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
147      REAL run_off_lic_0(klon), y_run_off_lic_0(klon)      ! ffonte----Flux thermique utilise pour fondre la neige
148        ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
149      REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)      !           hauteur de neige, en kg/m2/s
150      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)  
151    
152      REAL rain_f(klon), snow_f(klon)      REAL flux_o(klon), flux_g(klon)
153      REAL fder(klon)      !IM "slab" ocean
154        ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')
155        ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')
156    
157      REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)      REAL tslab(klon)
158      REAL rugos(klon, nbsrf)      ! tslab-in/output-R temperature du slab ocean (en Kelvin)
159      ! la nouvelle repartition des surfaces sortie de l'interface      ! uniqmnt pour slab
     REAL pctsrf_new(klon, nbsrf)  
160    
161      REAL zcoefh(klon, klev)      ! Local:
     REAL zu1(klon)  
     REAL zv1(klon)  
162    
163      !$$$ PB ajout pour soil      REAL y_flux_o(klon), y_flux_g(klon)
164      LOGICAL, INTENT (IN) :: soil_model      real ytslab(klon)
165      !IM ajout seuils cdrm, cdrh      REAL y_fqcalving(klon), y_ffonte(klon)
166      REAL cdmmax, cdhmax      real y_run_off_lic_0(klon)
167    
168      REAL ksta, ksta_ter      REAL rugmer(klon)
     LOGICAL ok_kzmin  
169    
     REAL ftsoil(klon, nsoilmx, nbsrf)  
170      REAL ytsoil(klon, nsoilmx)      REAL ytsoil(klon, nsoilmx)
     REAL qsol(klon)  
   
     EXTERNAL clqh, clvent, coefkz, calbeta, cltrac  
171    
172      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
173      REAL yalb(klon)      REAL yalb(klon)
     REAL yalblw(klon)  
174      REAL yu1(klon), yv1(klon)      REAL yu1(klon), yv1(klon)
175      REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)      ! on rajoute en output yu1 et yv1 qui sont les vents dans
176      REAL yrain_f(klon), ysnow_f(klon)      ! la premiere couche
177      REAL ysollw(klon), ysolsw(klon), ysollwdown(klon)      REAL ysnow(klon), yqsurf(klon), yagesno(klon)
178      REAL yfder(klon), ytaux(klon), ytauy(klon)  
179        real yqsol(klon)
180        ! column-density of water in soil, in kg m-2
181    
182        REAL yrain_f(klon)
183        ! liquid water mass flux (kg/m2/s), positive down
184    
185        REAL ysnow_f(klon)
186        ! solid water mass flux (kg/m2/s), positive down
187    
188        REAL ysollw(klon), ysolsw(klon)
189        REAL yfder(klon)
190      REAL yrugm(klon), yrads(klon), yrugoro(klon)      REAL yrugm(klon), yrads(klon), yrugoro(klon)
191    
192      REAL yfluxlat(klon)      REAL yfluxlat(klon)
# Line 199  contains Line 197  contains
197      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)
198      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)
199      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
200      REAL ycoefh(klon, klev), ycoefm(klon, klev)      REAL coefh(klon, klev), coefm(klon, klev)
201      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
202      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
203      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)
204    
     LOGICAL ok_nonloc  
     PARAMETER (ok_nonloc=.FALSE.)  
205      REAL ycoefm0(klon, klev), ycoefh0(klon, klev)      REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
206    
     !IM 081204 hcl_Anne ? BEG  
207      REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)      REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
208      REAL ykmm(klon, klev+1), ykmn(klon, klev+1)      REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
209      REAL ykmq(klon, klev+1)      REAL ykmq(klon, klev+1)
210      REAL yq2(klon, klev+1), q2(klon, klev+1, nbsrf)      REAL yq2(klon, klev+1)
211      REAL q2diag(klon, klev+1)      REAL q2diag(klon, klev+1)
     !IM 081204 hcl_Anne ? END  
212    
213      REAL u1lay(klon), v1lay(klon)      REAL u1lay(klon), v1lay(klon)
214      REAL delp(klon, klev)      REAL delp(klon, klev)
# Line 223  contains Line 217  contains
217      INTEGER ni(klon), knon, j      INTEGER ni(klon), knon, j
218    
219      REAL pctsrf_pot(klon, nbsrf)      REAL pctsrf_pot(klon, nbsrf)
220      ! "pourcentage potentiel" pour tenir compte des éventuelles      ! "pourcentage potentiel" pour tenir compte des \'eventuelles
221      ! apparitions ou disparitions de la glace de mer      ! apparitions ou disparitions de la glace de mer
222    
223      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
224    
     ! 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)  
   
225      REAL yt2m(klon), yq2m(klon), yu10m(klon)      REAL yt2m(klon), yq2m(klon), yu10m(klon)
226      REAL yustar(klon)      REAL yustar(klon)
227      ! -- LOOP      ! -- LOOP
# Line 257  contains Line 231  contains
231      ! -- LOOP      ! -- LOOP
232    
233      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)  
234      REAL ypblh(klon)      REAL ypblh(klon)
235      REAL ylcl(klon)      REAL ylcl(klon)
236      REAL ycapcl(klon)      REAL ycapcl(klon)
# Line 279  contains Line 241  contains
241      REAL ytrmb1(klon)      REAL ytrmb1(klon)
242      REAL ytrmb2(klon)      REAL ytrmb2(klon)
243      REAL ytrmb3(klon)      REAL ytrmb3(klon)
     REAL y_cd_h(klon), y_cd_m(klon)  
244      REAL uzon(klon), vmer(klon)      REAL uzon(klon), vmer(klon)
245      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
246      REAL psfce(klon), patm(klon)      REAL psfce(klon), patm(klon)
# Line 291  contains Line 252  contains
252      LOGICAL zxli      LOGICAL zxli
253      PARAMETER (zxli=.FALSE.)      PARAMETER (zxli=.FALSE.)
254    
     REAL zt, zqs, zdelta, zcor  
     REAL t_coup  
     PARAMETER (t_coup=273.15)  
   
     CHARACTER (len=20) :: modname = 'clmain'  
   
255      !------------------------------------------------------------      !------------------------------------------------------------
256    
257      ytherm = 0.      ytherm = 0.
258    
     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  
   
259      DO k = 1, klev ! epaisseur de couche      DO k = 1, klev ! epaisseur de couche
260         DO i = 1, klon         DO i = 1, klon
261            delp(i, k) = paprs(i, k) - paprs(i, k+1)            delp(i, k) = paprs(i, k) - paprs(i, k+1)
# Line 354  contains Line 280  contains
280      yts = 0.      yts = 0.
281      ysnow = 0.      ysnow = 0.
282      yqsurf = 0.      yqsurf = 0.
     yalb = 0.  
     yalblw = 0.  
283      yrain_f = 0.      yrain_f = 0.
284      ysnow_f = 0.      ysnow_f = 0.
285      yfder = 0.      yfder = 0.
     ytaux = 0.  
     ytauy = 0.  
286      ysolsw = 0.      ysolsw = 0.
287      ysollw = 0.      ysollw = 0.
     ysollwdown = 0.  
288      yrugos = 0.      yrugos = 0.
289      yu1 = 0.      yu1 = 0.
290      yv1 = 0.      yv1 = 0.
# Line 378  contains Line 299  contains
299      pctsrf_new = 0.      pctsrf_new = 0.
300      y_flux_u = 0.      y_flux_u = 0.
301      y_flux_v = 0.      y_flux_v = 0.
     !$$ PB  
302      y_dflux_t = 0.      y_dflux_t = 0.
303      y_dflux_q = 0.      y_dflux_q = 0.
304      ytsoil = 999999.      ytsoil = 999999.
305      yrugoro = 0.      yrugoro = 0.
     ! -- LOOP  
306      yu10mx = 0.      yu10mx = 0.
307      yu10my = 0.      yu10my = 0.
308      ywindsp = 0.      ywindsp = 0.
     ! -- LOOP  
309      d_ts = 0.      d_ts = 0.
     !§§§ PB  
310      yfluxlat = 0.      yfluxlat = 0.
311      flux_t = 0.      flux_t = 0.
312      flux_q = 0.      flux_q = 0.
# Line 399  contains Line 316  contains
316      d_q = 0.      d_q = 0.
317      d_u = 0.      d_u = 0.
318      d_v = 0.      d_v = 0.
319      zcoefh = 0.      ycoefh = 0.
   
     ! Boucler sur toutes les sous-fractions du sol:  
320    
321      ! Initialisation des "pourcentages potentiels". On considère ici qu'on      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
322      ! peut avoir potentiellement de la glace sur tout le domaine océanique      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
323      ! (à affiner)      ! (\`a affiner)
324    
325      pctsrf_pot = pctsrf      pctsrf_pot = pctsrf
326      pctsrf_pot(:, is_oce) = 1. - zmasq      pctsrf_pot(:, is_oce) = 1. - zmasq
327      pctsrf_pot(:, is_sic) = 1. - zmasq      pctsrf_pot(:, is_sic) = 1. - zmasq
328    
329      DO nsrf = 1, nbsrf      ! Boucler sur toutes les sous-fractions du sol:
330         ! chercher les indices:  
331        loop_surface: DO nsrf = 1, nbsrf
332           ! Chercher les indices :
333         ni = 0         ni = 0
334         knon = 0         knon = 0
335         DO i = 1, klon         DO i = 1, klon
336            ! Pour déterminer le domaine à traiter, on utilise les surfaces            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
337            ! "potentielles"            ! "potentielles"
338            IF (pctsrf_pot(i, nsrf) > epsfra) THEN            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
339               knon = knon + 1               knon = knon + 1
# Line 424  contains Line 341  contains
341            END IF            END IF
342         END DO         END DO
343    
344         ! 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  
345            DO j = 1, knon            DO j = 1, knon
346               i = ni(j)               i = ni(j)
347               ytsoil(j, k) = ftsoil(i, k, nsrf)               ypct(j) = pctsrf(i, nsrf)
348            END DO               yts(j) = ts(i, nsrf)
349         END DO               ytslab(i) = tslab(i)
350         DO k = 1, klev               ysnow(j) = snow(i, nsrf)
351            DO j = 1, knon               yqsurf(j) = qsurf(i, nsrf)
352               i = ni(j)               yalb(j) = falbe(i, nsrf)
353               ypaprs(j, k) = paprs(i, k)               yrain_f(j) = rain_fall(i)
354               ypplay(j, k) = pplay(i, k)               ysnow_f(j) = snow_f(i)
355               ydelp(j, k) = delp(i, k)               yagesno(j) = agesno(i, nsrf)
356               yu(j, k) = u(i, k)               yfder(j) = fder(i)
357               yv(j, k) = v(i, k)               ysolsw(j) = solsw(i, nsrf)
358               yt(j, k) = t(i, k)               ysollw(j) = sollw(i, nsrf)
359               yq(j, k) = q(i, k)               yrugos(j) = rugos(i, nsrf)
360            END DO               yrugoro(j) = rugoro(i)
361         END DO               yu1(j) = u1lay(i)
362                 yv1(j) = v1lay(i)
363                 yrads(j) = ysolsw(j) + ysollw(j)
364                 ypaprs(j, klev+1) = paprs(i, klev+1)
365                 y_run_off_lic_0(j) = run_off_lic_0(i)
366                 yu10mx(j) = u10m(i, nsrf)
367                 yu10my(j) = v10m(i, nsrf)
368                 ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))
369              END DO
370    
371              ! For continent, copy soil water content
372              IF (nsrf == is_ter) THEN
373                 yqsol(:knon) = qsol(ni(:knon))
374              ELSE
375                 yqsol = 0.
376              END IF
377    
378         ! calculer Cdrag et les coefficients d'echange            DO k = 1, nsoilmx
379         CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,&               DO j = 1, knon
380              yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)                  i = ni(j)
381         !IM 081204 BEG                  ytsoil(j, k) = ftsoil(i, k, nsrf)
        !CR test  
        IF (iflag_pbl==1) THEN  
           !IM 081204 END  
           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))  
382               END DO               END DO
383            END DO            END DO
        END IF  
   
        !IM cf JLD : on seuille ycoefm et ycoefh  
        IF (nsrf==is_oce) THEN  
           DO j = 1, knon  
              !           ycoefm(j, 1)=min(ycoefm(j, 1), 1.1E-3)  
              ycoefm(j, 1) = min(ycoefm(j, 1), cdmmax)  
              !           ycoefh(j, 1)=min(ycoefh(j, 1), 1.1E-3)  
              ycoefh(j, 1) = min(ycoefh(j, 1), cdhmax)  
           END DO  
        END IF  
   
        !IM: 261103  
        IF (ok_kzmin) THEN  
           !IM cf FH: 201103 BEG  
           !   Calcul d'une diffusion minimale pour les conditions tres stables.  
           CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, ycoefm, &  
                ycoefm0, ycoefh0)  
384    
           IF (1==1) THEN  
              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  
           !IM cf FH: 201103 END  
           !IM: 261103  
        END IF !ok_kzmin  
   
        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  
385            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  
386               DO j = 1, knon               DO j = 1, knon
387                  i = ni(j)                  i = ni(j)
388                  yq2(j, k) = q2(i, k, nsrf)                  ypaprs(j, k) = paprs(i, k)
389                    ypplay(j, k) = pplay(i, k)
390                    ydelp(j, k) = delp(i, k)
391                    yu(j, k) = u(i, k)
392                    yv(j, k) = v(i, k)
393                    yt(j, k) = t(i, k)
394                    yq(j, k) = q(i, k)
395               END DO               END DO
396            END DO            END DO
397    
398            !   Bug introduit volontairement pour converger avec les resultats            ! calculer Cdrag et les coefficients d'echange
399            !  du papier sur les thermiques.            CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &
400            IF (1==1) THEN                 yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))
401               y_cd_m(1:knon) = ycoefm(1:knon, 1)            IF (iflag_pbl == 1) THEN
402               y_cd_h(1:knon) = ycoefh(1:knon, 1)               CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
403            ELSE               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
404               y_cd_h(1:knon) = ycoefm(1:knon, 1)               coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
              y_cd_m(1:knon) = ycoefh(1:knon, 1)  
405            END IF            END IF
           CALL ustarhb(knon, yu, yv, y_cd_m, yustar)  
406    
407            IF (prt_level>9) THEN            ! on met un seuil pour coefm et coefh
408               PRINT *, 'USTAR = ', yustar            IF (nsrf == is_oce) THEN
409                 coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
410                 coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
411            END IF            END IF
412    
413            !   iflag_pbl peut etre utilise comme longuer de melange            IF (ok_kzmin) THEN
414                 ! Calcul d'une diffusion minimale pour les conditions tres stables
415                 CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
416                      coefm(:knon, 1), ycoefm0, ycoefh0)
417                 coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
418                 coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
419              END IF
420    
421            IF (iflag_pbl>=11) THEN            IF (iflag_pbl >= 3) THEN
422               CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &               ! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et
423                    yu, yv, yteta, y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, &               ! Fr\'ed\'eric Hourdin
424                    iflag_pbl)               yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
425            ELSE                    + ypplay(:knon, 1))) &
426               CALL yamada4(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, &                    * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
427                    yv, yteta, y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)               DO k = 2, klev
428                    yzlay(1:knon, k) = yzlay(1:knon, k-1) &
429                         + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
430                         / ypaprs(1:knon, k) &
431                         * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
432                 END DO
433                 DO k = 1, klev
434                    yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
435                         / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
436                 END DO
437                 yzlev(1:knon, 1) = 0.
438                 yzlev(:knon, klev+1) = 2. * yzlay(:knon, klev) &
439                      - yzlay(:knon, klev - 1)
440                 DO k = 2, klev
441                    yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
442                 END DO
443                 DO k = 1, klev + 1
444                    DO j = 1, knon
445                       i = ni(j)
446                       yq2(j, k) = q2(i, k, nsrf)
447                    END DO
448                 END DO
449    
450                 CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
451                 IF (prt_level > 9) PRINT *, 'USTAR = ', yustar
452    
453                 ! iflag_pbl peut \^etre utilis\'e comme longueur de m\'elange
454    
455                 IF (iflag_pbl >= 11) THEN
456                    CALL vdif_kcay(knon, dtime, rg, ypaprs, yzlev, yzlay, yu, yv, &
457                         yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, yustar, &
458                         iflag_pbl)
459                 ELSE
460                    CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
461                         coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
462                 END IF
463    
464                 coefm(:knon, 2:) = ykmm(:knon, 2:klev)
465                 coefh(:knon, 2:) = ykmn(:knon, 2:klev)
466            END IF            END IF
467    
468            ycoefm(1:knon, 1) = y_cd_m(1:knon)            ! calculer la diffusion des vitesses "u" et "v"
469            ycoefh(1:knon, 1) = y_cd_h(1:knon)            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &
470            ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)                 ypplay, ydelp, y_d_u, y_flux_u)
471            ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
472         END IF                 ypplay, ydelp, y_d_v, y_flux_v)
473    
474         ! calculer la diffusion des vitesses "u" et "v"            ! calculer la diffusion de "q" et de "h"
475         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &            CALL clqh(dtime, itap, jour, debut, rlat, knon, nsrf, ni(:knon), &
476              ydelp, y_d_u, y_flux_u)                 pctsrf, ytsoil, yqsol, rmu0, co2_ppm, yrugos, yrugoro, yu1, &
477         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &                 yv1, coefh(:knon, :), yt, yq, yts, ypaprs, ypplay, ydelp, &
478              ydelp, y_d_v, y_flux_v)                 yrads, yalb(:knon), ysnow, yqsurf, yrain_f, ysnow_f, yfder, &
479                   ysolsw, yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, &
480         ! pour le couplage                 y_d_ts(:knon), yz0_new, y_flux_t, y_flux_q, y_dflux_t, &
481         ytaux = y_flux_u(:, 1)                 y_dflux_q, y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, &
482         ytauy = y_flux_v(:, 1)                 y_flux_g)
483    
484         ! calculer la diffusion de "q" et de "h"            ! calculer la longueur de rugosite sur ocean
485         CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,&            yrugm = 0.
486              cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil,&            IF (nsrf == is_oce) THEN
487              yqsol, ok_veget, ocean, npas, nexca, rmu0, co2_ppm, yrugos,&               DO j = 1, knon
488              yrugoro, yu1, yv1, ycoefh, yt, yq, yts, ypaprs, ypplay,&                  yrugm(j) = 0.018*coefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
489              ydelp, yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, &                       0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))
490              yfder, ytaux, ytauy, ywindsp, ysollw, ysollwdown, ysolsw,&                  yrugm(j) = max(1.5E-05, yrugm(j))
491              yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts,&               END DO
492              yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,&            END IF
             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  
493            DO j = 1, knon            DO j = 1, knon
494               yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &               y_dflux_t(j) = y_dflux_t(j)*ypct(j)
495                    0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2))               y_dflux_q(j) = y_dflux_q(j)*ypct(j)
496               yrugm(j) = max(1.5E-05, yrugm(j))               yu1(j) = yu1(j)*ypct(j)
497            END DO               yv1(j) = yv1(j)*ypct(j)
498         END IF            END DO
        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  
499    
500         DO k = 1, klev            DO k = 1, klev
501            DO j = 1, knon               DO j = 1, knon
502               i = ni(j)                  i = ni(j)
503               ycoefh(j, k) = ycoefh(j, k)*ypct(j)                  coefh(j, k) = coefh(j, k)*ypct(j)
504               ycoefm(j, k) = ycoefm(j, k)*ypct(j)                  coefm(j, k) = coefm(j, k)*ypct(j)
505               y_d_t(j, k) = y_d_t(j, k)*ypct(j)                  y_d_t(j, k) = y_d_t(j, k)*ypct(j)
506               y_d_q(j, k) = y_d_q(j, k)*ypct(j)                  y_d_q(j, k) = y_d_q(j, k)*ypct(j)
507               !§§§ PB                  flux_t(i, k, nsrf) = y_flux_t(j, k)
508               flux_t(i, k, nsrf) = y_flux_t(j, k)                  flux_q(i, k, nsrf) = y_flux_q(j, k)
509               flux_q(i, k, nsrf) = y_flux_q(j, k)                  flux_u(i, k, nsrf) = y_flux_u(j, k)
510               flux_u(i, k, nsrf) = y_flux_u(j, k)                  flux_v(i, k, nsrf) = y_flux_v(j, k)
511               flux_v(i, k, nsrf) = y_flux_v(j, k)                  y_d_u(j, k) = y_d_u(j, k)*ypct(j)
512               !$$$ PB        y_flux_t(j, k) = y_flux_t(j, k) * ypct(j)                  y_d_v(j, k) = y_d_v(j, k)*ypct(j)
513               !$$$ PB        y_flux_q(j, k) = y_flux_q(j, k) * ypct(j)               END DO
              y_d_u(j, k) = y_d_u(j, k)*ypct(j)  
              y_d_v(j, k) = y_d_v(j, k)*ypct(j)  
              !$$$ PB        y_flux_u(j, k) = y_flux_u(j, k) * ypct(j)  
              !$$$ PB        y_flux_v(j, k) = y_flux_v(j, k) * ypct(j)  
514            END DO            END DO
        END DO  
515    
516         evap(:, nsrf) = -flux_q(:, 1, nsrf)            evap(:, nsrf) = -flux_q(:, 1, nsrf)
517    
518         albe(:, nsrf) = 0.            falbe(:, nsrf) = 0.
519         alblw(:, nsrf) = 0.            snow(:, nsrf) = 0.
520         snow(:, nsrf) = 0.            qsurf(:, nsrf) = 0.
521         qsurf(:, nsrf) = 0.            rugos(:, nsrf) = 0.
522         rugos(:, nsrf) = 0.            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)  
           !$$$ pb         rugmer(i) = yrugm(j)  
           IF (nsrf==is_oce) THEN  
              rugmer(i) = yrugm(j)  
              rugos(i, nsrf) = yrugm(j)  
           END IF  
           !IM cf JLD ??  
           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  
523            DO j = 1, knon            DO j = 1, knon
524               i = ni(j)               i = ni(j)
525               qsol(i) = yqsol(j)               d_ts(i, nsrf) = y_d_ts(j)
526                 falbe(i, nsrf) = yalb(j)
527                 snow(i, nsrf) = ysnow(j)
528                 qsurf(i, nsrf) = yqsurf(j)
529                 rugos(i, nsrf) = yz0_new(j)
530                 fluxlat(i, nsrf) = yfluxlat(j)
531                 IF (nsrf == is_oce) THEN
532                    rugmer(i) = yrugm(j)
533                    rugos(i, nsrf) = yrugm(j)
534                 END IF
535                 agesno(i, nsrf) = yagesno(j)
536                 fqcalving(i, nsrf) = y_fqcalving(j)
537                 ffonte(i, nsrf) = y_ffonte(j)
538                 cdragh(i) = cdragh(i) + coefh(j, 1)
539                 cdragm(i) = cdragm(i) + coefm(j, 1)
540                 dflux_t(i) = dflux_t(i) + y_dflux_t(j)
541                 dflux_q(i) = dflux_q(i) + y_dflux_q(j)
542                 zu1(i) = zu1(i) + yu1(j)
543                 zv1(i) = zv1(i) + yv1(j)
544              END DO
545              IF (nsrf == is_ter) THEN
546                 qsol(ni(:knon)) = yqsol(:knon)
547              else IF (nsrf == is_lic) THEN
548                 DO j = 1, knon
549                    i = ni(j)
550                    run_off_lic_0(i) = y_run_off_lic_0(j)
551                 END DO
552              END IF
553    
554              ftsoil(:, :, nsrf) = 0.
555              DO k = 1, nsoilmx
556                 DO j = 1, knon
557                    i = ni(j)
558                    ftsoil(i, k, nsrf) = ytsoil(j, k)
559                 END DO
560            END DO            END DO
561         END IF  
        IF (nsrf==is_lic) THEN  
562            DO j = 1, knon            DO j = 1, knon
563               i = ni(j)               i = ni(j)
564               run_off_lic_0(i) = y_run_off_lic_0(j)               DO k = 1, klev
565                    d_t(i, k) = d_t(i, k) + y_d_t(j, k)
566                    d_q(i, k) = d_q(i, k) + y_d_q(j, k)
567                    d_u(i, k) = d_u(i, k) + y_d_u(j, k)
568                    d_v(i, k) = d_v(i, k) + y_d_v(j, k)
569                    ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
570                 END DO
571            END DO            END DO
572         END IF  
573         !$$$ PB ajout pour soil            ! diagnostic t, q a 2m et u, v a 10m
574         ftsoil(:, :, nsrf) = 0.  
        DO k = 1, nsoilmx  
575            DO j = 1, knon            DO j = 1, knon
576               i = ni(j)               i = ni(j)
577               ftsoil(i, k, nsrf) = ytsoil(j, k)               uzon(j) = yu(j, 1) + y_d_u(j, 1)
578            END DO               vmer(j) = yv(j, 1) + y_d_v(j, 1)
579         END DO               tair1(j) = yt(j, 1) + y_d_t(j, 1)
580                 qair1(j) = yq(j, 1) + y_d_q(j, 1)
581                 zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
582                      1)))*(ypaprs(j, 1)-ypplay(j, 1))
583                 tairsol(j) = yts(j) + y_d_ts(j)
584                 rugo1(j) = yrugos(j)
585                 IF (nsrf == is_oce) THEN
586                    rugo1(j) = rugos(i, nsrf)
587                 END IF
588                 psfce(j) = ypaprs(j, 1)
589                 patm(j) = ypplay(j, 1)
590    
591         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)  
              !$$$ PB        flux_t(i, k) = flux_t(i, k) + y_flux_t(j, k)  
              !$$$         flux_q(i, k) = flux_q(i, k) + y_flux_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)  
              !$$$  PB       flux_u(i, k) = flux_u(i, k) + y_flux_u(j, k)  
              !$$$         flux_v(i, k) = flux_v(i, k) + y_flux_v(j, k)  
              zcoefh(i, k) = zcoefh(i, k) + ycoefh(j, k)  
592            END DO            END DO
        END DO  
593    
594         !cc diagnostic t, q a 2m et u, v a 10m            CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &
595                   zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &
596                   yt10m, yq10m, yu10m, yustar)
597    
598         DO j = 1, knon            DO j = 1, knon
599            i = ni(j)               i = ni(j)
600            uzon(j) = yu(j, 1) + y_d_u(j, 1)               t2m(i, nsrf) = yt2m(j)
601            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  
602    
603         CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &               ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
604              tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &               u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
605              yu10m, yustar)               v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)
        !IM 081204 END  
   
        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)  
606    
607         END DO            END DO
608    
609         DO i = 1, knon            CALL hbtm(knon, ypaprs, ypplay, yt2m, yq2m, yustar, &
610            y_cd_h(i) = ycoefh(i, 1)                 y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &
611            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  
612    
        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  
613            DO j = 1, knon            DO j = 1, knon
              ! on projette sur la grille globale  
614               i = ni(j)               i = ni(j)
615               IF (pctsrf_new(i, is_oce)>epsfra) THEN               pblh(i, nsrf) = ypblh(j)
616                  flux_o(i) = y_flux_o(j)               plcl(i, nsrf) = ylcl(j)
617               ELSE               capcl(i, nsrf) = ycapcl(j)
618                  flux_o(i) = 0.               oliqcl(i, nsrf) = yoliqcl(j)
619               END IF               cteicl(i, nsrf) = ycteicl(j)
620                 pblt(i, nsrf) = ypblt(j)
621                 therm(i, nsrf) = ytherm(j)
622                 trmb1(i, nsrf) = ytrmb1(j)
623                 trmb2(i, nsrf) = ytrmb2(j)
624                 trmb3(i, nsrf) = ytrmb3(j)
625            END DO            END DO
        END IF  
626    
        IF (nsrf==is_sic) THEN  
627            DO j = 1, knon            DO j = 1, knon
628               i = ni(j)               DO k = 1, klev + 1
629               ! On pondère lorsque l'on fait le bilan au sol :                  i = ni(j)
630               ! flux_g(i) = y_flux_g(j)*ypct(j)                  q2(i, k, nsrf) = yq2(j, k)
631               IF (pctsrf_new(i, is_sic)>epsfra) THEN               END DO
                 flux_g(i) = y_flux_g(j)  
              ELSE  
                 flux_g(i) = 0.  
              END IF  
632            END DO            END DO
633              !IM "slab" ocean
634              IF (nsrf == is_oce) THEN
635                 DO j = 1, knon
636                    ! on projette sur la grille globale
637                    i = ni(j)
638                    IF (pctsrf_new(i, is_oce)>epsfra) THEN
639                       flux_o(i) = y_flux_o(j)
640                    ELSE
641                       flux_o(i) = 0.
642                    END IF
643                 END DO
644              END IF
645    
646              IF (nsrf == is_sic) THEN
647                 DO j = 1, knon
648                    i = ni(j)
649                    ! On pond\`ere lorsque l'on fait le bilan au sol :
650                    IF (pctsrf_new(i, is_sic)>epsfra) THEN
651                       flux_g(i) = y_flux_g(j)
652                    ELSE
653                       flux_g(i) = 0.
654                    END IF
655                 END DO
656    
        END IF  
        !nsrf.EQ.is_sic                                              
        IF (ocean=='slab  ') THEN  
           IF (nsrf==is_oce) THEN  
              tslab(1:klon) = ytslab(1:klon)  
              seaice(1:klon) = y_seaice(1:klon)  
              !nsrf                                                        
657            END IF            END IF
658            !OCEAN                                                               end IF if_knon
659         END IF      END DO loop_surface
     END DO  
660    
661      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
     ! A rajouter: conservation de l'albedo  
662    
663      rugos(:, is_oce) = rugmer      rugos(:, is_oce) = rugmer
664      pctsrf = pctsrf_new      pctsrf = pctsrf_new

Legend:
Removed from v.40  
changed lines
  Added in v.155

  ViewVC Help
Powered by ViewVC 1.1.21