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

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

  ViewVC Help
Powered by ViewVC 1.1.21