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

Legend:
Removed from v.13  
changed lines
  Added in v.145

  ViewVC Help
Powered by ViewVC 1.1.21