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

Diff of /trunk/Sources/phylmd/clmain.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

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

  ViewVC Help
Powered by ViewVC 1.1.21