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

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

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

revision 38 by guez, Thu Jan 6 17:52:19 2011 UTC revision 62 by guez, Thu Jul 26 14:37:37 2012 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, zcoefh, 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 : "zcoefh",
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, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)
73        REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
74        REAL, INTENT(IN):: rlon(klon)
75        REAL, INTENT(IN):: rlat(klon) ! latitude en degrés
76        REAL cufi(klon), cvfi(klon)
77      ! cufi-----input-R- resolution des mailles en x (m)      ! cufi-----input-R- resolution des mailles en x (m)
78      ! cvfi-----input-R- resolution des mailles en y (m)      ! cvfi-----input-R- resolution des mailles en y (m)
79        REAL d_t(klon, klev), d_q(klon, klev)
80      ! d_t------output-R- le changement pour "t"      ! d_t------output-R- le changement pour "t"
81      ! d_q------output-R- le changement pour "q"      ! d_q------output-R- le changement pour "q"
82      ! d_u------output-R- le changement pour "u"  
83      ! d_v------output-R- le changement pour "v"      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
84      ! d_ts-----output-R- le changement pour "ts"      ! changement pour "u" et "v"
85    
86        REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)
87      ! 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)
88      !                    (orientation positive vers le bas)      !                    (orientation positive vers le bas)
89      ! 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)
90      ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal      REAL dflux_t(klon), dflux_q(klon)
     ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal  
91      ! dflux_t derive du flux sensible      ! dflux_t derive du flux sensible
92      ! dflux_q derive du flux latent      ! dflux_q derive du flux latent
93      !IM "slab" ocean      !IM "slab" ocean
94        REAL flux_o(klon), flux_g(klon)
95        !IM "slab" ocean
96      ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')      ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')
97      ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')      ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')
98        REAL y_flux_o(klon), y_flux_g(klon)
99        REAL tslab(klon), ytslab(klon)
100      ! tslab-in/output-R temperature du slab ocean (en Kelvin)      ! tslab-in/output-R temperature du slab ocean (en Kelvin)
101      ! uniqmnt pour slab      ! uniqmnt pour slab
102        REAL seaice(klon), y_seaice(klon)
103      ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')      ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')
104      !cc      REAL y_fqcalving(klon), y_ffonte(klon)
105        REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
106      ! ffonte----Flux thermique utilise pour fondre la neige      ! ffonte----Flux thermique utilise pour fondre la neige
107      ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la      ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
108      !           hauteur de neige, en kg/m2/s      !           hauteur de neige, en kg/m2/s
     ! on rajoute en output yu1 et yv1 qui sont les vents dans  
     ! la premiere couche  
     ! ces 4 variables sont maintenant traites dans phytrac  
     ! itr--------input-I- nombre de traceurs  
     ! tr---------input-R- q. de traceurs  
     ! flux_surf--input-R- flux de traceurs a la surface  
     ! d_tr-------output-R tendance de traceurs  
     !IM cf. AM : PBL  
     ! trmb1-------deep_cape  
     ! trmb2--------inhibition  
     ! trmb3-------Point Omega  
     ! Cape(klon)-------Cape du thermique  
     ! 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)  
     REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)  
109      REAL run_off_lic_0(klon), y_run_off_lic_0(klon)      REAL run_off_lic_0(klon), y_run_off_lic_0(klon)
110    
111      REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)      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
113        ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
114      REAL rugmer(klon), agesno(klon, nbsrf)      REAL rugmer(klon), agesno(klon, nbsrf)
115      REAL, INTENT (IN) :: rugoro(klon)      REAL, INTENT(IN):: rugoro(klon)
116      REAL cdragh(klon), cdragm(klon)      REAL, INTENT(out):: cdragh(klon), cdragm(klon)
     ! jour de l'annee en cours                  
     INTEGER jour  
     REAL rmu0(klon) ! cosinus de l'angle solaire zenithal      
117      ! taux CO2 atmosphere                          ! taux CO2 atmosphere                    
118      REAL co2_ppm      REAL co2_ppm
119      LOGICAL, INTENT (IN) :: debut      LOGICAL, INTENT(IN):: debut
120      LOGICAL, INTENT (IN) :: lafin      LOGICAL, INTENT(IN):: lafin
121      LOGICAL ok_veget      LOGICAL ok_veget
122      CHARACTER (len=*), INTENT (IN) :: ocean      CHARACTER(len=*), INTENT(IN):: ocean
123      INTEGER npas, nexca      INTEGER npas, nexca
124    
     REAL pctsrf(klon, nbsrf)  
125      REAL ts(klon, nbsrf)      REAL ts(klon, nbsrf)
126        ! ts-------input-R- temperature du sol (en Kelvin)
127      REAL d_ts(klon, nbsrf)      REAL d_ts(klon, nbsrf)
128        ! d_ts-----output-R- le changement pour "ts"
129      REAL snow(klon, nbsrf)      REAL snow(klon, nbsrf)
130      REAL qsurf(klon, nbsrf)      REAL qsurf(klon, nbsrf)
131      REAL evap(klon, nbsrf)      REAL evap(klon, nbsrf)
# Line 155  contains Line 134  contains
134    
135      REAL fluxlat(klon, nbsrf)      REAL fluxlat(klon, nbsrf)
136    
137      REAL rain_f(klon), snow_f(klon)      REAL, intent(in):: rain_fall(klon), snow_f(klon)
138      REAL fder(klon)      REAL fder(klon)
139    
140      REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)      REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)
141      REAL rugos(klon, nbsrf)      REAL rugos(klon, nbsrf)
142      ! la nouvelle repartition des surfaces sortie de l'interface      ! rugos----input-R- longeur de rugosite (en m)
     REAL pctsrf_new(klon, nbsrf)  
143    
144      REAL zcoefh(klon, klev)      REAL zcoefh(klon, klev)
145      REAL zu1(klon)      REAL zu1(klon)
146      REAL zv1(klon)      REAL zv1(klon)
147    
148      !$$$ PB ajout pour soil      !$$$ PB ajout pour soil
149      LOGICAL, INTENT (IN) :: soil_model      LOGICAL, INTENT(IN):: soil_model
150      !IM ajout seuils cdrm, cdrh      !IM ajout seuils cdrm, cdrh
151      REAL cdmmax, cdhmax      REAL cdmmax, cdhmax
152    
# Line 179  contains Line 157  contains
157      REAL ytsoil(klon, nsoilmx)      REAL ytsoil(klon, nsoilmx)
158      REAL qsol(klon)      REAL qsol(klon)
159    
     EXTERNAL clqh, clvent, coefkz, calbeta, cltrac  
   
160      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
161      REAL yalb(klon)      REAL yalb(klon)
162      REAL yalblw(klon)      REAL yalblw(klon)
163      REAL yu1(klon), yv1(klon)      REAL yu1(klon), yv1(klon)
164        ! on rajoute en output yu1 et yv1 qui sont les vents dans
165        ! la premiere couche
166      REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)      REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)
167      REAL yrain_f(klon), ysnow_f(klon)      REAL yrain_f(klon), ysnow_f(klon)
168      REAL ysollw(klon), ysolsw(klon), ysollwdown(klon)      REAL ysollw(klon), ysolsw(klon), ysollwdown(klon)
# Line 199  contains Line 177  contains
177      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)
178      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)
179      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
180      REAL ycoefh(klon, klev), ycoefm(klon, klev)      REAL coefh(klon, klev), coefm(klon, klev)
181      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
182      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
183      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 186  contains
186      PARAMETER (ok_nonloc=.FALSE.)      PARAMETER (ok_nonloc=.FALSE.)
187      REAL ycoefm0(klon, klev), ycoefh0(klon, klev)      REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
188    
     !IM 081204 hcl_Anne ? BEG  
189      REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)      REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
190      REAL ykmm(klon, klev+1), ykmn(klon, klev+1)      REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
191      REAL ykmq(klon, klev+1)      REAL ykmq(klon, klev+1)
192      REAL yq2(klon, klev+1), q2(klon, klev+1, nbsrf)      REAL yq2(klon, klev+1), q2(klon, klev+1, nbsrf)
193      REAL q2diag(klon, klev+1)      REAL q2diag(klon, klev+1)
     !IM 081204 hcl_Anne ? END  
194    
195      REAL u1lay(klon), v1lay(klon)      REAL u1lay(klon), v1lay(klon)
196      REAL delp(klon, klev)      REAL delp(klon, klev)
197      INTEGER i, k, nsrf      INTEGER i, k, nsrf
198    
199      INTEGER ni(klon), knon, j      INTEGER ni(klon), knon, j
200      ! Introduction d'une variable "pourcentage potentiel" pour tenir compte  
     ! des eventuelles apparitions et/ou disparitions de la glace de mer  
201      REAL pctsrf_pot(klon, nbsrf)      REAL pctsrf_pot(klon, nbsrf)
202        ! "pourcentage potentiel" pour tenir compte des éventuelles
203        ! apparitions ou disparitions de la glace de mer
204    
205      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
206    
207      ! maf pour sorties IOISPL en cas de debugagage      ! maf pour sorties IOISPL en cas de debugagage
208    
209      CHARACTER (80) cldebug      CHARACTER(80) cldebug
210      SAVE cldebug      SAVE cldebug
211      CHARACTER (8) cl_surf(nbsrf)      CHARACTER(8) cl_surf(nbsrf)
212      SAVE cl_surf      SAVE cl_surf
213      INTEGER nhoridbg, nidbg      INTEGER nhoridbg, nidbg
214      SAVE nhoridbg, nidbg      SAVE nhoridbg, nidbg
# Line 242  contains Line 219  contains
219      LOGICAL first_appel      LOGICAL first_appel
220      SAVE first_appel      SAVE first_appel
221      DATA first_appel/ .TRUE./      DATA first_appel/ .TRUE./
222      LOGICAL :: debugindex = .FALSE.      LOGICAL:: debugindex = .FALSE.
223      INTEGER idayref      INTEGER idayref
224      REAL t2m(klon, nbsrf), q2m(klon, nbsrf)      REAL t2m(klon, nbsrf), q2m(klon, nbsrf)
225      REAL u10m(klon, nbsrf), v10m(klon, nbsrf)      REAL u10m(klon, nbsrf), v10m(klon, nbsrf)
# Line 259  contains Line 236  contains
236      !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds      !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds
237      ! physiq ce qui permet de sortir les grdeurs par sous surface)      ! physiq ce qui permet de sortir les grdeurs par sous surface)
238      REAL pblh(klon, nbsrf)      REAL pblh(klon, nbsrf)
239        ! pblh------- HCL
240      REAL plcl(klon, nbsrf)      REAL plcl(klon, nbsrf)
241      REAL capcl(klon, nbsrf)      REAL capcl(klon, nbsrf)
242      REAL oliqcl(klon, nbsrf)      REAL oliqcl(klon, nbsrf)
243      REAL cteicl(klon, nbsrf)      REAL cteicl(klon, nbsrf)
244      REAL pblt(klon, nbsrf)      REAL pblt(klon, nbsrf)
245        ! pblT------- T au nveau HCL
246      REAL therm(klon, nbsrf)      REAL therm(klon, nbsrf)
247      REAL trmb1(klon, nbsrf)      REAL trmb1(klon, nbsrf)
248        ! trmb1-------deep_cape
249      REAL trmb2(klon, nbsrf)      REAL trmb2(klon, nbsrf)
250        ! trmb2--------inhibition
251      REAL trmb3(klon, nbsrf)      REAL trmb3(klon, nbsrf)
252        ! trmb3-------Point Omega
253      REAL ypblh(klon)      REAL ypblh(klon)
254      REAL ylcl(klon)      REAL ylcl(klon)
255      REAL ycapcl(klon)      REAL ycapcl(klon)
# Line 278  contains Line 260  contains
260      REAL ytrmb1(klon)      REAL ytrmb1(klon)
261      REAL ytrmb2(klon)      REAL ytrmb2(klon)
262      REAL ytrmb3(klon)      REAL ytrmb3(klon)
     REAL y_cd_h(klon), y_cd_m(klon)  
263      REAL uzon(klon), vmer(klon)      REAL uzon(klon), vmer(klon)
264      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
265      REAL psfce(klon), patm(klon)      REAL psfce(klon), patm(klon)
# Line 294  contains Line 275  contains
275      REAL t_coup      REAL t_coup
276      PARAMETER (t_coup=273.15)      PARAMETER (t_coup=273.15)
277    
278      CHARACTER (len=20) :: modname = 'clmain'      CHARACTER(len=20):: modname = 'clmain'
279    
280      !------------------------------------------------------------      !------------------------------------------------------------
281    
     ! initialisation Anne  
282      ytherm = 0.      ytherm = 0.
283    
284      IF (debugindex .AND. first_appel) THEN      IF (debugindex .AND. first_appel) THEN
# Line 307  contains Line 287  contains
287         ! initialisation sorties netcdf         ! initialisation sorties netcdf
288    
289         idayref = day_ini         idayref = day_ini
290         CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)         CALL ymds2ju(annee_ref, 1, idayref, 0., zjulian)
291         CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlon, zx_lon)         CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlon, zx_lon)
292         DO i = 1, iim         DO i = 1, iim
293            zx_lon(i, 1) = rlon(i+1)            zx_lon(i, 1) = rlon(i+1)
# Line 342  contains Line 322  contains
322         v1lay(i) = v(i, 1)*zx_alf1 + v(i, 2)*zx_alf2         v1lay(i) = v(i, 1)*zx_alf1 + v(i, 2)*zx_alf2
323      END DO      END DO
324    
325      ! initialisation:      ! Initialization:
326        rugmer = 0.
327      DO i = 1, klon      cdragh = 0.
328         rugmer(i) = 0.0      cdragm = 0.
329         cdragh(i) = 0.0      dflux_t = 0.
330         cdragm(i) = 0.0      dflux_q = 0.
331         dflux_t(i) = 0.0      zu1 = 0.
332         dflux_q(i) = 0.0      zv1 = 0.
333         zu1(i) = 0.0      ypct = 0.
334         zv1(i) = 0.0      yts = 0.
335      END DO      ysnow = 0.
336      ypct = 0.0      yqsurf = 0.
337      yts = 0.0      yalb = 0.
338      ysnow = 0.0      yalblw = 0.
339      yqsurf = 0.0      yrain_f = 0.
340      yalb = 0.0      ysnow_f = 0.
341      yalblw = 0.0      yfder = 0.
342      yrain_f = 0.0      ytaux = 0.
343      ysnow_f = 0.0      ytauy = 0.
344      yfder = 0.0      ysolsw = 0.
345      ytaux = 0.0      ysollw = 0.
346      ytauy = 0.0      ysollwdown = 0.
347      ysolsw = 0.0      yrugos = 0.
348      ysollw = 0.0      yu1 = 0.
349      ysollwdown = 0.0      yv1 = 0.
350      yrugos = 0.0      yrads = 0.
351      yu1 = 0.0      ypaprs = 0.
352      yv1 = 0.0      ypplay = 0.
353      yrads = 0.0      ydelp = 0.
354      ypaprs = 0.0      yu = 0.
355      ypplay = 0.0      yv = 0.
356      ydelp = 0.0      yt = 0.
357      yu = 0.0      yq = 0.
358      yv = 0.0      pctsrf_new = 0.
359      yt = 0.0      y_flux_u = 0.
360      yq = 0.0      y_flux_v = 0.
     pctsrf_new = 0.0  
     y_flux_u = 0.0  
     y_flux_v = 0.0  
361      !$$ PB      !$$ PB
362      y_dflux_t = 0.0      y_dflux_t = 0.
363      y_dflux_q = 0.0      y_dflux_q = 0.
364      ytsoil = 999999.      ytsoil = 999999.
365      yrugoro = 0.      yrugoro = 0.
366      ! -- LOOP      ! -- LOOP
367      yu10mx = 0.0      yu10mx = 0.
368      yu10my = 0.0      yu10my = 0.
369      ywindsp = 0.0      ywindsp = 0.
370      ! -- LOOP      ! -- LOOP
371      DO nsrf = 1, nbsrf      d_ts = 0.
        DO i = 1, klon  
           d_ts(i, nsrf) = 0.0  
        END DO  
     END DO  
372      !§§§ PB      !§§§ PB
373      yfluxlat = 0.      yfluxlat = 0.
374      flux_t = 0.      flux_t = 0.
375      flux_q = 0.      flux_q = 0.
376      flux_u = 0.      flux_u = 0.
377      flux_v = 0.      flux_v = 0.
378      DO k = 1, klev      d_t = 0.
379         DO i = 1, klon      d_q = 0.
380            d_t(i, k) = 0.0      d_u = 0.
381            d_q(i, k) = 0.0      d_v = 0.
382            d_u(i, k) = 0.0      zcoefh = 0.
           d_v(i, k) = 0.0  
           zcoefh(i, k) = 0.0  
        END DO  
     END DO  
383    
384      ! Boucler sur toutes les sous-fractions du sol:      ! Boucler sur toutes les sous-fractions du sol:
385    
# Line 422  contains Line 391  contains
391      pctsrf_pot(:, is_oce) = 1. - zmasq      pctsrf_pot(:, is_oce) = 1. - zmasq
392      pctsrf_pot(:, is_sic) = 1. - zmasq      pctsrf_pot(:, is_sic) = 1. - zmasq
393    
394      DO nsrf = 1, nbsrf      loop_surface: DO nsrf = 1, nbsrf
395         ! chercher les indices:         ! Chercher les indices :
396         ni = 0         ni = 0
397         knon = 0         knon = 0
398         DO i = 1, klon         DO i = 1, klon
399            ! pour determiner le domaine a traiter on utilise les surfaces            ! Pour déterminer le domaine à traiter, on utilise les surfaces
400            ! "potentielles"            ! "potentielles"
401            IF (pctsrf_pot(i, nsrf) > epsfra) THEN            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
402               knon = knon + 1               knon = knon + 1
# Line 447  contains Line 416  contains
416            CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)            CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)
417         END IF         END IF
418    
419         IF (knon==0) CYCLE         if_knon: IF (knon /= 0) then
   
        DO j = 1, knon  
           i = ni(j)  
           ypct(j) = pctsrf(i, nsrf)  
           yts(j) = ts(i, nsrf)  
           ytslab(i) = tslab(i)  
           ysnow(j) = snow(i, nsrf)  
           yqsurf(j) = qsurf(i, nsrf)  
           yalb(j) = albe(i, nsrf)  
           yalblw(j) = alblw(i, nsrf)  
           yrain_f(j) = rain_f(i)  
           ysnow_f(j) = snow_f(i)  
           yagesno(j) = agesno(i, nsrf)  
           yfder(j) = fder(i)  
           ytaux(j) = flux_u(i, 1, nsrf)  
           ytauy(j) = flux_v(i, 1, nsrf)  
           ysolsw(j) = solsw(i, nsrf)  
           ysollw(j) = sollw(i, nsrf)  
           ysollwdown(j) = sollwdown(i)  
           yrugos(j) = rugos(i, nsrf)  
           yrugoro(j) = rugoro(i)  
           yu1(j) = u1lay(i)  
           yv1(j) = v1lay(i)  
           yrads(j) = ysolsw(j) + ysollw(j)  
           ypaprs(j, klev+1) = paprs(i, klev+1)  
           y_run_off_lic_0(j) = run_off_lic_0(i)  
           yu10mx(j) = u10m(i, nsrf)  
           yu10my(j) = v10m(i, nsrf)  
           ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))  
        END DO  
   
        !     IF bucket model for continent, copy soil water content  
        IF (nsrf==is_ter .AND. .NOT. ok_veget) THEN  
           DO j = 1, knon  
              i = ni(j)  
              yqsol(j) = qsol(i)  
           END DO  
        ELSE  
           yqsol = 0.  
        END IF  
        !$$$ PB ajour pour soil  
        DO k = 1, nsoilmx  
           DO j = 1, knon  
              i = ni(j)  
              ytsoil(j, k) = ftsoil(i, k, nsrf)  
           END DO  
        END DO  
        DO k = 1, klev  
420            DO j = 1, knon            DO j = 1, knon
421               i = ni(j)               i = ni(j)
422               ypaprs(j, k) = paprs(i, k)               ypct(j) = pctsrf(i, nsrf)
423               ypplay(j, k) = pplay(i, k)               yts(j) = ts(i, nsrf)
424               ydelp(j, k) = delp(i, k)               ytslab(i) = tslab(i)
425               yu(j, k) = u(i, k)               ysnow(j) = snow(i, nsrf)
426               yv(j, k) = v(i, k)               yqsurf(j) = qsurf(i, nsrf)
427               yt(j, k) = t(i, k)               yalb(j) = albe(i, nsrf)
428               yq(j, k) = q(i, k)               yalblw(j) = alblw(i, nsrf)
429                 yrain_f(j) = rain_fall(i)
430                 ysnow_f(j) = snow_f(i)
431                 yagesno(j) = agesno(i, nsrf)
432                 yfder(j) = fder(i)
433                 ytaux(j) = flux_u(i, 1, nsrf)
434                 ytauy(j) = flux_v(i, 1, nsrf)
435                 ysolsw(j) = solsw(i, nsrf)
436                 ysollw(j) = sollw(i, nsrf)
437                 ysollwdown(j) = sollwdown(i)
438                 yrugos(j) = rugos(i, nsrf)
439                 yrugoro(j) = rugoro(i)
440                 yu1(j) = u1lay(i)
441                 yv1(j) = v1lay(i)
442                 yrads(j) = ysolsw(j) + ysollw(j)
443                 ypaprs(j, klev+1) = paprs(i, klev+1)
444                 y_run_off_lic_0(j) = run_off_lic_0(i)
445                 yu10mx(j) = u10m(i, nsrf)
446                 yu10my(j) = v10m(i, nsrf)
447                 ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))
448            END DO            END DO
        END DO  
449    
450         ! calculer Cdrag et les coefficients d'echange            ! IF bucket model for continent, copy soil water content
451         CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,&            IF (nsrf == is_ter .AND. .NOT. ok_veget) THEN
452              yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)               DO j = 1, knon
453         !IM 081204 BEG                  i = ni(j)
454         !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))  
455               END DO               END DO
456            END DO            ELSE
457         END IF               yqsol = 0.
458              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)  
459    
460            IF (1==1) THEN            DO k = 1, nsoilmx
461               DO k = 1, klev               DO j = 1, knon
462                  DO i = 1, knon                  i = ni(j)
463                     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  
464               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  
465            END DO            END DO
466    
467            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  
468               DO j = 1, knon               DO j = 1, knon
469                  i = ni(j)                  i = ni(j)
470                  yq2(j, k) = q2(i, k, nsrf)                  ypaprs(j, k) = paprs(i, k)
471                    ypplay(j, k) = pplay(i, k)
472                    ydelp(j, k) = delp(i, k)
473                    yu(j, k) = u(i, k)
474                    yv(j, k) = v(i, k)
475                    yt(j, k) = t(i, k)
476                    yq(j, k) = q(i, k)
477               END DO               END DO
478            END DO            END DO
479    
480            !   Bug introduit volontairement pour converger avec les resultats            ! calculer Cdrag et les coefficients d'echange
481            !  du papier sur les thermiques.            CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &
482            IF (1==1) THEN                 yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))
483               y_cd_m(1:knon) = ycoefm(1:knon, 1)            IF (iflag_pbl == 1) THEN
484               y_cd_h(1:knon) = ycoefh(1:knon, 1)               CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
485            ELSE               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
486               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)  
487            END IF            END IF
           CALL ustarhb(knon, yu, yv, y_cd_m, yustar)  
488    
489            IF (prt_level>9) THEN            ! on seuille coefm et coefh
490               PRINT *, 'USTAR = ', yustar            IF (nsrf == is_oce) THEN
491                 coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
492                 coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
493            END IF            END IF
494    
495            !   iflag_pbl peut etre utilise comme longuer de melange            IF (ok_kzmin) THEN
496                 ! Calcul d'une diffusion minimale pour les conditions tres stables
497                 CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
498                      coefm(:, 1), ycoefm0, ycoefh0)
499                 coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
500                 coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
501               END IF
502    
503              IF (iflag_pbl >= 3) THEN
504                 ! MELLOR ET YAMADA adapté à Mars, Richard Fournier et
505                 ! Frédéric Hourdin
506                 yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
507                      + ypplay(:knon, 1))) &
508                      * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
509                 DO k = 2, klev
510                    yzlay(1:knon, k) = yzlay(1:knon, k-1) &
511                         + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
512                         / ypaprs(1:knon, k) &
513                         * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
514                 END DO
515                 DO k = 1, klev
516                    yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
517                         / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
518                 END DO
519                 yzlev(1:knon, 1) = 0.
520                 yzlev(:knon, klev+1) = 2. * yzlay(:knon, klev) &
521                      - yzlay(:knon, klev - 1)
522                 DO k = 2, klev
523                    yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
524                 END DO
525                 DO k = 1, klev + 1
526                    DO j = 1, knon
527                       i = ni(j)
528                       yq2(j, k) = q2(i, k, nsrf)
529                    END DO
530                 END DO
531    
532                 CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
533    
534            IF (iflag_pbl>=11) THEN               IF (prt_level > 9) THEN
535               CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &                  PRINT *, 'USTAR = ', yustar
536                    yu, yv, yteta, y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, &               END IF
                   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  
537    
538            ycoefm(1:knon, 1) = y_cd_m(1:knon)               ! iflag_pbl peut être utilisé comme longueur de mélange
539            ycoefh(1:knon, 1) = y_cd_h(1:knon)  
540            ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)               IF (iflag_pbl >= 11) THEN
541            ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)                  CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &
542         END IF                       yu, yv, yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, &
543                         yustar, iflag_pbl)
544                 ELSE
545                    CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
546                         coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
547                 END IF
548    
549                 coefm(:knon, 2:) = ykmm(:knon, 2:klev)
550                 coefh(:knon, 2:) = ykmn(:knon, 2:klev)
551              END IF
552    
553         ! calculer la diffusion des vitesses "u" et "v"            ! calculer la diffusion des vitesses "u" et "v"
554         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &            CALL clvent(knon, dtime, yu1, yv1, coefm, yt, yu, ypaprs, ypplay, &
555              ydelp, y_d_u, y_flux_u)                 ydelp, y_d_u, y_flux_u)
556         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &            CALL clvent(knon, dtime, yu1, yv1, coefm, yt, yv, ypaprs, ypplay, &
557              ydelp, y_d_v, y_flux_v)                 ydelp, y_d_v, y_flux_v)
558    
559         ! pour le couplage            ! pour le couplage
560         ytaux = y_flux_u(:, 1)            ytaux = y_flux_u(:, 1)
561         ytauy = y_flux_v(:, 1)            ytauy = y_flux_v(:, 1)
562    
563         ! FH modif sur le cdrag temperature            ! calculer la diffusion de "q" et de "h"
564         !$$$PB : déplace dans clcdrag            CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat, &
565         !$$$      do i=1, knon                 cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil, &
566         !$$$         ycoefh(i, 1)=ycoefm(i, 1)*0.8                 yqsol, ok_veget, ocean, npas, nexca, rmu0, co2_ppm, yrugos, &
567         !$$$      enddo                 yrugoro, yu1, yv1, coefh, yt, yq, yts, ypaprs, ypplay, &
568                   ydelp, yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, &
569         ! calculer la diffusion de "q" et de "h"                 yfder, ytaux, ytauy, ywindsp, ysollw, ysollwdown, ysolsw, &
570         CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,&                 yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts, &
571              cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil,&                 yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q, &
572              yqsol, ok_veget, ocean, npas, nexca, rmu0, co2_ppm, yrugos,&                 y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g, &
573              yrugoro, yu1, yv1, ycoefh, yt, yq, yts, ypaprs, ypplay,&                 ytslab, y_seaice)
574              ydelp, yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, &  
575              yfder, ytaux, ytauy, ywindsp, ysollw, ysollwdown, ysolsw,&            ! calculer la longueur de rugosite sur ocean
576              yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts,&            yrugm = 0.
577              yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,&            IF (nsrf == is_oce) THEN
578              y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g,&               DO j = 1, knon
579              ytslab, y_seaice)                  yrugm(j) = 0.018*coefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
580                         0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))
581         ! calculer la longueur de rugosite sur ocean                  yrugm(j) = max(1.5E-05, yrugm(j))
582         yrugm = 0.               END DO
583         IF (nsrf==is_oce) THEN            END IF
584            DO j = 1, knon            DO j = 1, knon
585               yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &               y_dflux_t(j) = y_dflux_t(j)*ypct(j)
586                    0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2))               y_dflux_q(j) = y_dflux_q(j)*ypct(j)
587               yrugm(j) = max(1.5E-05, yrugm(j))               yu1(j) = yu1(j)*ypct(j)
588                 yv1(j) = yv1(j)*ypct(j)
589            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  
590    
591         DO k = 1, klev            DO k = 1, klev
592            DO j = 1, knon               DO j = 1, knon
593               i = ni(j)                  i = ni(j)
594               ycoefh(j, k) = ycoefh(j, k)*ypct(j)                  coefh(j, k) = coefh(j, k)*ypct(j)
595               ycoefm(j, k) = ycoefm(j, k)*ypct(j)                  coefm(j, k) = coefm(j, k)*ypct(j)
596               y_d_t(j, k) = y_d_t(j, k)*ypct(j)                  y_d_t(j, k) = y_d_t(j, k)*ypct(j)
597               y_d_q(j, k) = y_d_q(j, k)*ypct(j)                  y_d_q(j, k) = y_d_q(j, k)*ypct(j)
598               !§§§ PB                  flux_t(i, k, nsrf) = y_flux_t(j, k)
599               flux_t(i, k, nsrf) = y_flux_t(j, k)                  flux_q(i, k, nsrf) = y_flux_q(j, k)
600               flux_q(i, k, nsrf) = y_flux_q(j, k)                  flux_u(i, k, nsrf) = y_flux_u(j, k)
601               flux_u(i, k, nsrf) = y_flux_u(j, k)                  flux_v(i, k, nsrf) = y_flux_v(j, k)
602               flux_v(i, k, nsrf) = y_flux_v(j, k)                  y_d_u(j, k) = y_d_u(j, k)*ypct(j)
603               !$$$ 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)
604               !$$$ 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)  
605            END DO            END DO
        END DO  
606    
607         evap(:, nsrf) = -flux_q(:, 1, nsrf)            evap(:, nsrf) = -flux_q(:, 1, nsrf)
608    
609         albe(:, nsrf) = 0.            albe(:, nsrf) = 0.
610         alblw(:, nsrf) = 0.            alblw(:, nsrf) = 0.
611         snow(:, nsrf) = 0.            snow(:, nsrf) = 0.
612         qsurf(:, nsrf) = 0.            qsurf(:, nsrf) = 0.
613         rugos(:, nsrf) = 0.            rugos(:, nsrf) = 0.
614         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  
615            DO j = 1, knon            DO j = 1, knon
616               i = ni(j)               i = ni(j)
617               qsol(i) = yqsol(j)               d_ts(i, nsrf) = y_d_ts(j)
618                 albe(i, nsrf) = yalb(j)
619                 alblw(i, nsrf) = yalblw(j)
620                 snow(i, nsrf) = ysnow(j)
621                 qsurf(i, nsrf) = yqsurf(j)
622                 rugos(i, nsrf) = yz0_new(j)
623                 fluxlat(i, nsrf) = yfluxlat(j)
624                 IF (nsrf == is_oce) THEN
625                    rugmer(i) = yrugm(j)
626                    rugos(i, nsrf) = yrugm(j)
627                 END IF
628                 agesno(i, nsrf) = yagesno(j)
629                 fqcalving(i, nsrf) = y_fqcalving(j)
630                 ffonte(i, nsrf) = y_ffonte(j)
631                 cdragh(i) = cdragh(i) + coefh(j, 1)
632                 cdragm(i) = cdragm(i) + coefm(j, 1)
633                 dflux_t(i) = dflux_t(i) + y_dflux_t(j)
634                 dflux_q(i) = dflux_q(i) + y_dflux_q(j)
635                 zu1(i) = zu1(i) + yu1(j)
636                 zv1(i) = zv1(i) + yv1(j)
637            END DO            END DO
638         END IF            IF (nsrf == is_ter) THEN
639         IF (nsrf==is_lic) THEN               DO j = 1, knon
640                    i = ni(j)
641                    qsol(i) = yqsol(j)
642                 END DO
643              END IF
644              IF (nsrf == is_lic) THEN
645                 DO j = 1, knon
646                    i = ni(j)
647                    run_off_lic_0(i) = y_run_off_lic_0(j)
648                 END DO
649              END IF
650              !$$$ PB ajout pour soil
651              ftsoil(:, :, nsrf) = 0.
652              DO k = 1, nsoilmx
653                 DO j = 1, knon
654                    i = ni(j)
655                    ftsoil(i, k, nsrf) = ytsoil(j, k)
656                 END DO
657              END DO
658    
659            DO j = 1, knon            DO j = 1, knon
660               i = ni(j)               i = ni(j)
661               run_off_lic_0(i) = y_run_off_lic_0(j)               DO k = 1, klev
662                    d_t(i, k) = d_t(i, k) + y_d_t(j, k)
663                    d_q(i, k) = d_q(i, k) + y_d_q(j, k)
664                    d_u(i, k) = d_u(i, k) + y_d_u(j, k)
665                    d_v(i, k) = d_v(i, k) + y_d_v(j, k)
666                    zcoefh(i, k) = zcoefh(i, k) + coefh(j, k)
667                 END DO
668            END DO            END DO
669         END IF  
670         !$$$ PB ajout pour soil            !cc diagnostic t, q a 2m et u, v a 10m
671         ftsoil(:, :, nsrf) = 0.  
        DO k = 1, nsoilmx  
672            DO j = 1, knon            DO j = 1, knon
673               i = ni(j)               i = ni(j)
674               ftsoil(i, k, nsrf) = ytsoil(j, k)               uzon(j) = yu(j, 1) + y_d_u(j, 1)
675            END DO               vmer(j) = yv(j, 1) + y_d_v(j, 1)
676         END DO               tair1(j) = yt(j, 1) + y_d_t(j, 1)
677                 qair1(j) = yq(j, 1) + y_d_q(j, 1)
678                 zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
679                      1)))*(ypaprs(j, 1)-ypplay(j, 1))
680                 tairsol(j) = yts(j) + y_d_ts(j)
681                 rugo1(j) = yrugos(j)
682                 IF (nsrf == is_oce) THEN
683                    rugo1(j) = rugos(i, nsrf)
684                 END IF
685                 psfce(j) = ypaprs(j, 1)
686                 patm(j) = ypplay(j, 1)
687    
688         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)  
689            END DO            END DO
        END DO  
690    
691         !cc diagnostic t, q a 2m et u, v a 10m            CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &
692                   zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &
693                   yt10m, yq10m, yu10m, yustar)
694    
695         DO j = 1, knon            DO j = 1, knon
696            i = ni(j)               i = ni(j)
697            uzon(j) = yu(j, 1) + y_d_u(j, 1)               t2m(i, nsrf) = yt2m(j)
698            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  
699    
700         CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &               ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
701              tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &               u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
702              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)  
703    
704         END DO            END DO
705    
706         DO i = 1, knon            CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &
707            y_cd_h(i) = ycoefh(i, 1)                 y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &
708            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  
709    
        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  
710            DO j = 1, knon            DO j = 1, knon
              ! on projette sur la grille globale  
711               i = ni(j)               i = ni(j)
712               IF (pctsrf_new(i, is_oce)>epsfra) THEN               pblh(i, nsrf) = ypblh(j)
713                  flux_o(i) = y_flux_o(j)               plcl(i, nsrf) = ylcl(j)
714               ELSE               capcl(i, nsrf) = ycapcl(j)
715                  flux_o(i) = 0.               oliqcl(i, nsrf) = yoliqcl(j)
716               END IF               cteicl(i, nsrf) = ycteicl(j)
717                 pblt(i, nsrf) = ypblt(j)
718                 therm(i, nsrf) = ytherm(j)
719                 trmb1(i, nsrf) = ytrmb1(j)
720                 trmb2(i, nsrf) = ytrmb2(j)
721                 trmb3(i, nsrf) = ytrmb3(j)
722            END DO            END DO
        END IF  
723    
        IF (nsrf==is_sic) THEN  
724            DO j = 1, knon            DO j = 1, knon
725               i = ni(j)               DO k = 1, klev + 1
726               ! On pondère lorsque l'on fait le bilan au sol :                  i = ni(j)
727               ! flux_g(i) = y_flux_g(j)*ypct(j)                  q2(i, k, nsrf) = yq2(j, k)
728               IF (pctsrf_new(i, is_sic)>epsfra) THEN               END DO
                 flux_g(i) = y_flux_g(j)  
              ELSE  
                 flux_g(i) = 0.  
              END IF  
729            END DO            END DO
730              !IM "slab" ocean
731              IF (nsrf == is_oce) THEN
732                 DO j = 1, knon
733                    ! on projette sur la grille globale
734                    i = ni(j)
735                    IF (pctsrf_new(i, is_oce)>epsfra) THEN
736                       flux_o(i) = y_flux_o(j)
737                    ELSE
738                       flux_o(i) = 0.
739                    END IF
740                 END DO
741              END IF
742    
743              IF (nsrf == is_sic) THEN
744                 DO j = 1, knon
745                    i = ni(j)
746                    ! On pondère lorsque l'on fait le bilan au sol :
747                    IF (pctsrf_new(i, is_sic)>epsfra) THEN
748                       flux_g(i) = y_flux_g(j)
749                    ELSE
750                       flux_g(i) = 0.
751                    END IF
752                 END DO
753    
        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                                                        
754            END IF            END IF
755            !OCEAN                                                                  IF (ocean == 'slab  ') THEN
756         END IF               IF (nsrf == is_oce) THEN
757      END DO                  tslab(1:klon) = ytslab(1:klon)
758                    seaice(1:klon) = y_seaice(1:klon)
759                 END IF
760              END IF
761           end IF if_knon
762        END DO loop_surface
763    
764      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
     ! A rajouter: conservation de l'albedo  
765    
766      rugos(:, is_oce) = rugmer      rugos(:, is_oce) = rugmer
767      pctsrf = pctsrf_new      pctsrf = pctsrf_new

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

  ViewVC Help
Powered by ViewVC 1.1.21