/[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 10 by guez, Fri Apr 18 14:45:53 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, intent(in):: 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, intent(in):: 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           zx_alf1 = 1.0                  ypaprs(j, k) = paprs(i, k)
360           zx_alf2 = 1.0 - zx_alf1                  ypplay(j, k) = pplay(i, k)
361           u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2                  ydelp(j, k) = delp(i, k)
362           v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2                  yu(j, k) = u(i, k)
363        ENDDO                  yv(j, k) = v(i, k)
364  c                  yt(j, k) = t(i, k)
365  c initialisation:                  yq(j, k) = q(i, k)
366  c               END DO
367        DO i = 1, klon            END DO
368           rugmer(i) = 0.0  
369           cdragh(i) = 0.0            ! calculer Cdrag et les coefficients d'echange
370           cdragm(i) = 0.0            CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &
371           dflux_t(i) = 0.0                 yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))
372           dflux_q(i) = 0.0            IF (iflag_pbl == 1) THEN
373           zu1(i) = 0.0               CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
374           zv1(i) = 0.0               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
375        ENDDO               coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
376        ypct = 0.0            END IF
377        yts = 0.0  
378        ysnow = 0.0            ! on met un seuil pour coefm et coefh
379        yqsurf = 0.0            IF (nsrf == is_oce) THEN
380        yalb = 0.0               coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
381        yalblw = 0.0               coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
382        yrain_f = 0.0            END IF
383        ysnow_f = 0.0  
384        yfder = 0.0            IF (ok_kzmin) THEN
385        ytaux = 0.0               ! Calcul d'une diffusion minimale pour les conditions tres stables
386        ytauy = 0.0               CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
387        ysolsw = 0.0                    coefm(:knon, 1), ycoefm0, ycoefh0)
388        ysollw = 0.0               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
389        ysollwdown = 0.0               coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
390        yrugos = 0.0            END IF
391        yu1 = 0.0  
392        yv1 = 0.0            IF (iflag_pbl >= 3) THEN
393        yrads = 0.0               ! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et
394        ypaprs = 0.0               ! Fr\'ed\'eric Hourdin
395        ypplay = 0.0               yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
396        ydelp = 0.0                    + ypplay(:knon, 1))) &
397        yu = 0.0                    * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
398        yv = 0.0               DO k = 2, klev
399        yt = 0.0                  yzlay(1:knon, k) = yzlay(1:knon, k-1) &
400        yq = 0.0                       + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
401        pctsrf_new = 0.0                       / ypaprs(1:knon, k) &
402        y_flux_u = 0.0                       * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
403        y_flux_v = 0.0               END DO
404  C$$ PB               DO k = 1, klev
405        y_dflux_t = 0.0                  yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
406        y_dflux_q = 0.0                       / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
407        ytsoil = 999999.               END DO
408        yrugoro = 0.               yzlev(1:knon, 1) = 0.
409  c -- LOOP               yzlev(:knon, klev+1) = 2. * yzlay(:knon, klev) &
410        yu10mx = 0.0                    - yzlay(:knon, klev - 1)
411        yu10my = 0.0               DO k = 2, klev
412        ywindsp = 0.0                  yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
413  c -- LOOP               END DO
414        DO nsrf = 1, nbsrf               DO k = 1, klev + 1
415        DO i = 1, klon                  DO j = 1, knon
416           d_ts(i,nsrf) = 0.0                     i = ni(j)
417        ENDDO                     yq2(j, k) = q2(i, k, nsrf)
418        END DO                  END DO
419  C§§§ PB               END DO
420        yfluxlat=0.  
421        flux_t = 0.               CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
422        flux_q = 0.               IF (prt_level > 9) PRINT *, 'USTAR = ', yustar
423        flux_u = 0.  
424        flux_v = 0.               ! iflag_pbl peut \^etre utilis\'e comme longueur de m\'elange
425        DO k = 1, klev  
426        DO i = 1, klon               IF (iflag_pbl >= 11) THEN
427           d_t(i,k) = 0.0                  CALL vdif_kcay(knon, dtime, rg, ypaprs, yzlev, yzlay, yu, yv, &
428           d_q(i,k) = 0.0                       yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, yustar, &
429  c$$$         flux_t(i,k) = 0.0                       iflag_pbl)
430  c$$$         flux_q(i,k) = 0.0               ELSE
431           d_u(i,k) = 0.0                  CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
432           d_v(i,k) = 0.0                       coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
433  c$$$         flux_u(i,k) = 0.0               END IF
434  c$$$         flux_v(i,k) = 0.0  
435           zcoefh(i,k) = 0.0               coefm(:knon, 2:) = ykmm(:knon, 2:klev)
436        ENDDO               coefh(:knon, 2:) = ykmn(:knon, 2:klev)
437        ENDDO            END IF
438  cAA      IF (itr.GE.1) THEN  
439  cAA      DO it = 1, itr            ! calculer la diffusion des vitesses "u" et "v"
440  cAA      DO k = 1, klev            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &
441  cAA      DO i = 1, klon                 ypplay, ydelp, y_d_u, y_flux_u)
442  cAA         d_tr(i,k,it) = 0.0            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
443  cAA      ENDDO                 ypplay, ydelp, y_d_v, y_flux_v)
444  cAA      ENDDO  
445  cAA      ENDDO            ! calculer la diffusion de "q" et de "h"
446  cAA      ENDIF            CALL clqh(dtime, itap, jour, debut, rlat, knon, nsrf, ni(:knon), &
447                   pctsrf, ytsoil, yqsol, rmu0, yrugos, yrugoro, yu1, &
448  c                 yv1, coefh(:knon, :), yt, yq, yts, ypaprs, ypplay, ydelp, &
449  c Boucler sur toutes les sous-fractions du sol:                 yrads, yalb(:knon), ysnow, yqsurf, yrain_f, ysnow_f, yfder, &
450  c                 yfluxlat, pctsrf_new, yagesno(:knon), y_d_t, y_d_q, &
451  C Initialisation des "pourcentages potentiels". On considere ici qu'on                 y_d_ts(:knon), yz0_new, y_flux_t, y_flux_q, y_dflux_t, &
452  C peut avoir potentiellementdela glace sur tout le domaine oceanique                 y_dflux_q, y_fqcalving, y_ffonte, y_run_off_lic_0)
453  C (a affiner)  
454              ! calculer la longueur de rugosite sur ocean
455        pctsrf_pot = pctsrf            yrugm = 0.
456        pctsrf_pot(:,is_oce) = 1. - zmasq(:)            IF (nsrf == is_oce) THEN
457        pctsrf_pot(:,is_sic) = 1. - zmasq(:)               DO j = 1, knon
458                    yrugm(j) = 0.018*coefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
459        DO 99999 nsrf = 1, nbsrf                       0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))
460                    yrugm(j) = max(1.5E-05, yrugm(j))
461  c chercher les indices:               END DO
462        DO j = 1, klon            END IF
463           ni(j) = 0            DO j = 1, knon
464        ENDDO               y_dflux_t(j) = y_dflux_t(j)*ypct(j)
465        knon = 0               y_dflux_q(j) = y_dflux_q(j)*ypct(j)
466        DO i = 1, klon               yu1(j) = yu1(j)*ypct(j)
467                 yv1(j) = yv1(j)*ypct(j)
468  C pour determiner le domaine a traiter on utilise les surfaces "potentielles"            END DO
469  C    
470        IF (pctsrf_pot(i,nsrf).GT.epsfra) THEN            DO k = 1, klev
471           knon = knon + 1               DO j = 1, knon
472           ni(knon) = i                  i = ni(j)
473        ENDIF                  coefh(j, k) = coefh(j, k)*ypct(j)
474        ENDDO                  coefm(j, k) = coefm(j, k)*ypct(j)
475  c                  y_d_t(j, k) = y_d_t(j, k)*ypct(j)
476        if (check) THEN                  y_d_q(j, k) = y_d_q(j, k)*ypct(j)
477            write(*,*)'CLMAIN, nsrf, knon =',nsrf, knon                  flux_t(i, k, nsrf) = y_flux_t(j, k)
478  CC        call flush(6)                  flux_q(i, k, nsrf) = y_flux_q(j, k)
479        endif                  flux_u(i, k, nsrf) = y_flux_u(j, k)
480  c                  flux_v(i, k, nsrf) = y_flux_v(j, k)
481  c variables pour avoir une sortie IOIPSL des INDEX                  y_d_u(j, k) = y_d_u(j, k)*ypct(j)
482  c                  y_d_v(j, k) = y_d_v(j, k)*ypct(j)
483        IF (debugindex) THEN               END DO
484            tabindx(:)=0.            END DO
485  c          tabindx(1:knon)=(/FLOAT(i),i=1:knon/)  
486            DO i=1,knon            evap(:, nsrf) = -flux_q(:, 1, nsrf)
487              tabindx(1:knon)=FLOAT(i)  
488            END DO            falbe(:, nsrf) = 0.
489            debugtab(:,:)=0.            snow(:, nsrf) = 0.
490            ndexbg(:)=0            qsurf(:, nsrf) = 0.
491            CALL gath2cpl(tabindx,debugtab,klon,knon,iim,jjm,ni)            rugos(:, nsrf) = 0.
492            CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1)            fluxlat(:, nsrf) = 0.
493       $        ,ndexbg)            DO j = 1, knon
494        ENDIF               i = ni(j)
495        IF (knon.EQ.0) GOTO 99999               d_ts(i, nsrf) = y_d_ts(j)
496        DO j = 1, knon               falbe(i, nsrf) = yalb(j)
497        i = ni(j)               snow(i, nsrf) = ysnow(j)
498          ypct(j) = pctsrf(i,nsrf)               qsurf(i, nsrf) = yqsurf(j)
499          yts(j) = ts(i,nsrf)               rugos(i, nsrf) = yz0_new(j)
500  cIM "slab" ocean               fluxlat(i, nsrf) = yfluxlat(j)
501  c        PRINT *, 'tslab = ', i, tslab(i)               IF (nsrf == is_oce) THEN
502          ytslab(i) = tslab(i)                  rugmer(i) = yrugm(j)
503  c                  rugos(i, nsrf) = yrugm(j)
504          ysnow(j) = snow(i,nsrf)               END IF
505          yqsurf(j) = qsurf(i,nsrf)               agesno(i, nsrf) = yagesno(j)
506          yalb(j) = albe(i,nsrf)               fqcalving(i, nsrf) = y_fqcalving(j)
507          yalblw(j) = alblw(i,nsrf)               ffonte(i, nsrf) = y_ffonte(j)
508          yrain_f(j) = rain_f(i)               cdragh(i) = cdragh(i) + coefh(j, 1)
509          ysnow_f(j) = snow_f(i)               cdragm(i) = cdragm(i) + coefm(j, 1)
510          yagesno(j) = agesno(i,nsrf)               dflux_t(i) = dflux_t(i) + y_dflux_t(j)
511          yfder(j) = fder(i)               dflux_q(i) = dflux_q(i) + y_dflux_q(j)
512          ytaux(j) = flux_u(i,1,nsrf)               zu1(i) = zu1(i) + yu1(j)
513          ytauy(j) = flux_v(i,1,nsrf)               zv1(i) = zv1(i) + yv1(j)
514          ysolsw(j) = solsw(i,nsrf)            END DO
515          ysollw(j) = sollw(i,nsrf)            IF (nsrf == is_ter) THEN
516          ysollwdown(j) = sollwdown(i)               qsol(ni(:knon)) = yqsol(:knon)
517          yrugos(j) = rugos(i,nsrf)            else IF (nsrf == is_lic) THEN
518          yrugoro(j) = rugoro(i)               DO j = 1, knon
519          yu1(j) = u1lay(i)                  i = ni(j)
520          yv1(j) = v1lay(i)                  run_off_lic_0(i) = y_run_off_lic_0(j)
521          yrads(j) =  ysolsw(j)+ ysollw(j)               END DO
522          ypaprs(j,klev+1) = paprs(i,klev+1)            END IF
523          y_run_off_lic_0(j) = run_off_lic_0(i)  
524  c -- LOOP            ftsoil(:, :, nsrf) = 0.
525         yu10mx(j) = u10m(i,nsrf)            DO k = 1, nsoilmx
526         yu10my(j) = v10m(i,nsrf)               DO j = 1, knon
527         ywindsp(j) = SQRT(yu10mx(j)*yu10mx(j) + yu10my(j)*yu10my(j) )                  i = ni(j)
528  c -- LOOP                  ftsoil(i, k, nsrf) = ytsoil(j, k)
529        END DO               END DO
530  C            END DO
531  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.10  
changed lines
  Added in v.186

  ViewVC Help
Powered by ViewVC 1.1.21