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

Legend:
Removed from v.13  
changed lines
  Added in v.225

  ViewVC Help
Powered by ViewVC 1.1.21