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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21