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

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

  ViewVC Help
Powered by ViewVC 1.1.21