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

Diff of /trunk/phylmd/Interface_surf/pbl_surface.f

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

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

Legend:
Removed from v.3  
changed lines
  Added in v.186

  ViewVC Help
Powered by ViewVC 1.1.21