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

Diff of /trunk/phylmd/clmain.f

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21