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

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

  ViewVC Help
Powered by ViewVC 1.1.21