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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21