/[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 10 by guez, Fri Apr 18 14:45:53 2008 UTC trunk/Sources/phylmd/clmain.f revision 207 by guez, Thu Sep 1 10:30: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) ! 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 "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 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),rugoro(klon)      ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
145        REAL cdragh(klon), cdragm(klon)      !           hauteur de neige, en kg/m2/s
146        integer jour            ! jour de l'annee en cours      REAL run_off_lic_0(klon)
147        real rmu0(klon)         ! cosinus de l'angle solaire zenithal  
148        REAL co2_ppm            ! taux CO2 atmosphere      ! Local:
149        LOGICAL, intent(in):: debut  
150        logical, intent(in):: lafin      LOGICAL:: firstcal = .true.
151        logical ok_veget  
152        character*6 ocean      ! la nouvelle repartition des surfaces sortie de l'interface
153        integer npas, nexca      REAL, save:: pctsrf_new_oce(klon)
154  c      REAL, save:: pctsrf_new_sic(klon)
155        REAL pctsrf(klon,nbsrf)  
156        REAL ts(klon,nbsrf)      REAL y_fqcalving(klon), y_ffonte(klon)
157        REAL d_ts(klon,nbsrf)      real y_run_off_lic_0(klon)
158        REAL snow(klon,nbsrf)      REAL rugmer(klon)
159        REAL qsurf(klon,nbsrf)      REAL ytsoil(klon, nsoilmx)
160        REAL evap(klon,nbsrf)      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
161        REAL albe(klon,nbsrf)      REAL yalb(klon)
162        REAL alblw(klon,nbsrf)      REAL yu1(klon), yv1(klon)
163  c$$$ PB      ! on rajoute en output yu1 et yv1 qui sont les vents dans
164        REAL fluxlat(klon,nbsrf)      ! la premiere couche
165  C      REAL ysnow(klon), yqsurf(klon), yagesno(klon)
166        real rain_f(klon), snow_f(klon)  
167        REAL fder(klon)      real yqsol(klon)
168  cIM cf. JLD   REAL sollw(klon), solsw(klon), sollwdown(klon)      ! column-density of water in soil, in kg m-2
169        REAL sollw(klon,nbsrf), solsw(klon,nbsrf), sollwdown(klon)  
170        REAL rugos(klon,nbsrf)      REAL yrain_f(klon)
171  C la nouvelle repartition des surfaces sortie de l'interface      ! liquid water mass flux (kg/m2/s), positive down
172        REAL pctsrf_new(klon,nbsrf)  
173  cAA      REAL ysnow_f(klon)
174        REAL zcoefh(klon,klev)      ! solid water mass flux (kg/m2/s), positive down
175        REAL zu1(klon)  
176        REAL zv1(klon)      REAL yfder(klon)
177  cAA      REAL yrugm(klon), yrads(klon), yrugoro(klon)
178  c$$$ PB ajout pour soil  
179        LOGICAL soil_model      REAL yfluxlat(klon)
180  cIM ajout seuils cdrm, cdrh  
181        REAL cdmmax, cdhmax      REAL y_d_ts(klon)
182  cIM: 261103      REAL y_d_t(klon, klev), y_d_q(klon, klev)
183        REAL ksta, ksta_ter      REAL y_d_u(klon, klev), y_d_v(klon, klev)
184        LOGICAL ok_kzmin      REAL y_flux_t(klon), y_flux_q(klon)
185  cIM: 261103      REAL y_flux_u(klon), y_flux_v(klon)
186        REAL ftsoil(klon,nsoilmx,nbsrf)      REAL y_dflux_t(klon), y_dflux_q(klon)
187        REAL ytsoil(klon,nsoilmx)      REAL coefh(klon, klev), coefm(klon, klev)
188        REAL qsol(klon)      REAL yu(klon, klev), yv(klon, klev)
189  c======================================================================      REAL yt(klon, klev), yq(klon, klev)
190        EXTERNAL clqh, clvent, coefkz, calbeta, cltrac      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)
191  c======================================================================  
192        REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
193        REAL yalb(klon)  
194        REAL yalblw(klon)      REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
195        REAL yu1(klon), yv1(klon)      REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
196        real ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)      REAL ykmq(klon, klev+1)
197        real yrain_f(klon), ysnow_f(klon)      REAL yq2(klon, klev+1)
198        real ysollw(klon), ysolsw(klon), ysollwdown(klon)      REAL q2diag(klon, klev+1)
199        real yfder(klon), ytaux(klon), ytauy(klon)  
200        REAL yrugm(klon), yrads(klon),yrugoro(klon)      REAL u1lay(klon), v1lay(klon)
201  c$$$ PB      REAL delp(klon, klev)
202        REAL yfluxlat(klon)      INTEGER i, k, nsrf
203  C  
204        REAL y_d_ts(klon)      INTEGER ni(klon), knon, j
205        REAL y_d_t(klon, klev), y_d_q(klon, klev)  
206        REAL y_d_u(klon, klev), y_d_v(klon, klev)      REAL pctsrf_pot(klon, nbsrf)
207        REAL y_flux_t(klon,klev), y_flux_q(klon,klev)      ! "pourcentage potentiel" pour tenir compte des \'eventuelles
208        REAL y_flux_u(klon,klev), y_flux_v(klon,klev)      ! apparitions ou disparitions de la glace de mer
209        REAL y_dflux_t(klon), y_dflux_q(klon)  
210        REAL ycoefh(klon,klev), ycoefm(klon,klev)      REAL zx_alf1, zx_alf2 ! valeur ambiante par extrapolation
211        REAL yu(klon,klev), yv(klon,klev)  
212        REAL yt(klon,klev), yq(klon,klev)      REAL yt2m(klon), yq2m(klon), yu10m(klon)
213        REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)      REAL yustar(klon)
214  c  
215        LOGICAL ok_nonloc      REAL yt10m(klon), yq10m(klon)
216        PARAMETER (ok_nonloc=.FALSE.)      REAL ypblh(klon)
217        REAL ycoefm0(klon,klev), ycoefh0(klon,klev)      REAL ylcl(klon)
218        REAL ycapcl(klon)
219  cIM 081204 hcl_Anne ? BEG      REAL yoliqcl(klon)
220        real yzlay(klon,klev),yzlev(klon,klev+1),yteta(klon,klev)      REAL ycteicl(klon)
221        real ykmm(klon,klev+1),ykmn(klon,klev+1)      REAL ypblt(klon)
222        real ykmq(klon,klev+1)      REAL ytherm(klon)
223        real yq2(klon,klev+1),q2(klon,klev+1,nbsrf)      REAL ytrmb1(klon)
224        real q2diag(klon,klev+1)      REAL ytrmb2(klon)
225  cIM 081204   real yustar(klon),y_cd_m(klon),y_cd_h(klon)      REAL ytrmb3(klon)
226  cIM 081204 hcl_Anne ? END      REAL uzon(klon), vmer(klon)
227  c      REAL tair1(klon), qair1(klon), tairsol(klon)
228        REAL u1lay(klon), v1lay(klon)      REAL psfce(klon), patm(klon)
229        REAL delp(klon,klev)  
230        INTEGER i, k, nsrf      REAL qairsol(klon), zgeo1(klon)
231  cAA   INTEGER it      REAL rugo1(klon)
232        INTEGER ni(klon), knon, j  
233  c Introduction d'une variable "pourcentage potentiel" pour tenir compte      ! utiliser un jeu de fonctions simples              
234  c des eventuelles apparitions et/ou disparitions de la glace de mer      LOGICAL zxli
235        REAL pctsrf_pot(klon,nbsrf)      PARAMETER (zxli=.FALSE.)
236          
237  c======================================================================      !------------------------------------------------------------
238        REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.  
239  c======================================================================      ytherm = 0.
240  c  
241  c maf pour sorties IOISPL en cas de debugagage      DO k = 1, klev ! epaisseur de couche
242  c         DO i = 1, klon
243        CHARACTER*80 cldebug            delp(i, k) = paprs(i, k) - paprs(i, k+1)
244        SAVE cldebug         END DO
245        CHARACTER*8 cl_surf(nbsrf)      END DO
246        SAVE cl_surf      DO i = 1, klon ! vent de la premiere couche
247        INTEGER nhoridbg, nidbg         zx_alf1 = 1.0
248        SAVE nhoridbg, nidbg         zx_alf2 = 1.0 - zx_alf1
249        INTEGER ndexbg(iim*(jjm+1))         u1lay(i) = u(i, 1)*zx_alf1 + u(i, 2)*zx_alf2
250        REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1), zjulian         v1lay(i) = v(i, 1)*zx_alf1 + v(i, 2)*zx_alf2
251        REAL tabindx(klon)      END DO
252        REAL debugtab(iim,jjm+1)  
253        LOGICAL first_appel      ! Initialization:
254        SAVE first_appel      rugmer = 0.
255        DATA first_appel/.true./      cdragh = 0.
256        LOGICAL debugindex      cdragm = 0.
257        SAVE debugindex      dflux_t = 0.
258        DATA debugindex/.false./      dflux_q = 0.
259        integer idayref      zu1 = 0.
260        REAL t2m(klon,nbsrf), q2m(klon,nbsrf)      zv1 = 0.
261        REAL u10m(klon,nbsrf), v10m(klon,nbsrf)      ypct = 0.
262  c      yts = 0.
263        REAL yt2m(klon), yq2m(klon), yu10m(klon)      ysnow = 0.
264        REAL yustar(klon)      yqsurf = 0.
265  c -- LOOP      yrain_f = 0.
266         REAL yu10mx(klon)      ysnow_f = 0.
267         REAL yu10my(klon)      yfder = 0.
268         REAL ywindsp(klon)      yrugos = 0.
269  c -- LOOP      yu1 = 0.
270  c      yv1 = 0.
271        REAL yt10m(klon), yq10m(klon)      yrads = 0.
272  cIM cf. AM : pbl, hbtm2 (Comme les autres diagnostics on cumule ds physic ce qui      ypaprs = 0.
273  c   permet de sortir les grdeurs par sous surface)      ypplay = 0.
274        REAL pblh(klon,nbsrf)      ydelp = 0.
275        REAL plcl(klon,nbsrf)      yu = 0.
276        REAL capCL(klon,nbsrf)      yv = 0.
277        REAL oliqCL(klon,nbsrf)      yt = 0.
278        REAL cteiCL(klon,nbsrf)      yq = 0.
279        REAL pblT(klon,nbsrf)      y_dflux_t = 0.
280        REAL therm(klon,nbsrf)      y_dflux_q = 0.
281        REAL trmb1(klon,nbsrf)      ytsoil = 999999.
282        REAL trmb2(klon,nbsrf)      yrugoro = 0.
283        REAL trmb3(klon,nbsrf)      d_ts = 0.
284        REAL ypblh(klon)      yfluxlat = 0.
285        REAL ylcl(klon)      flux_t = 0.
286        REAL ycapCL(klon)      flux_q = 0.
287        REAL yoliqCL(klon)      flux_u = 0.
288        REAL ycteiCL(klon)      flux_v = 0.
289        REAL ypblT(klon)      d_t = 0.
290        REAL ytherm(klon)      d_q = 0.
291        REAL ytrmb1(klon)      d_u = 0.
292        REAL ytrmb2(klon)      d_v = 0.
293        REAL ytrmb3(klon)      ycoefh = 0.
294        REAL y_cd_h(klon), y_cd_m(klon)  
295  c     REAL ygamt(klon,2:klev) ! contre-gradient pour temperature      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
296  c     REAL ygamq(klon,2:klev) ! contre-gradient pour humidite      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
297        REAL uzon(klon), vmer(klon)      ! (\`a affiner)
298        REAL tair1(klon), qair1(klon), tairsol(klon)  
299        REAL psfce(klon), patm(klon)      pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
300  c      pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
301        REAL qairsol(klon), zgeo1(klon)      pctsrf_pot(:, is_oce) = 1. - zmasq
302        REAL rugo1(klon)      pctsrf_pot(:, is_sic) = 1. - zmasq
303  c  
304        LOGICAL zxli ! utiliser un jeu de fonctions simples      ! Tester si c'est le moment de lire le fichier:
305        PARAMETER (zxli=.FALSE.)      if (mod(itap - 1, lmt_pas) == 0) then
306  c         CALL interfoce_lim(jour, pctsrf_new_oce, pctsrf_new_sic)
307        REAL zt, zqs, zdelta, zcor      endif
308        REAL t_coup  
309        PARAMETER(t_coup=273.15)      ! Boucler sur toutes les sous-fractions du sol:
310  C  
311        character (len = 20) :: modname = 'clmain'      loop_surface: DO nsrf = 1, nbsrf
312        LOGICAL check         ! Chercher les indices :
313        PARAMETER (check=.false.)         ni = 0
314           knon = 0
315           DO i = 1, klon
316  c initialisation Anne            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
317        ytherm(:) = 0.            ! "potentielles"
318  C            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
319        if (check) THEN               knon = knon + 1
320            write(*,*) modname,'  klon=',klon               ni(knon) = i
321  CC        call flush(6)            END IF
322        endif         END DO
323        IF (debugindex .and. first_appel) THEN  
324            first_appel=.false.         if_knon: IF (knon /= 0) then
325  !            DO j = 1, knon
326  ! initialisation sorties netcdf               i = ni(j)
327  !               ypct(j) = pctsrf(i, nsrf)
328            idayref = day_ini               yts(j) = ftsol(i, nsrf)
329            CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)               ysnow(j) = snow(i, nsrf)
330            CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)               yqsurf(j) = qsurf(i, nsrf)
331            DO i = 1, iim               yalb(j) = falbe(i, nsrf)
332              zx_lon(i,1) = rlon(i+1)               yrain_f(j) = rain_fall(i)
333              zx_lon(i,jjm+1) = rlon(i+1)               ysnow_f(j) = snow_f(i)
334            ENDDO               yagesno(j) = agesno(i, nsrf)
335            CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)               yfder(j) = fder(i)
336            cldebug='sous_index'               yrugos(j) = rugos(i, nsrf)
337            CALL histbeg_totreg(cldebug, iim,zx_lon(:,1),jjm+1,               yrugoro(j) = rugoro(i)
338       $        zx_lat(1,:),1,iim,1,jjm               yu1(j) = u1lay(i)
339       $        +1, itau_phy,zjulian,dtime,nhoridbg,nidbg)               yv1(j) = v1lay(i)
340  ! no vertical axis               yrads(j) = solsw(i, nsrf) + sollw(i, nsrf)
341            cl_surf(1)='ter'               ypaprs(j, klev+1) = paprs(i, klev+1)
342            cl_surf(2)='lic'               y_run_off_lic_0(j) = run_off_lic_0(i)
343            cl_surf(3)='oce'            END DO
344            cl_surf(4)='sic'  
345            DO nsrf=1,nbsrf            ! For continent, copy soil water content
346              CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim,            IF (nsrf == is_ter) THEN
347       $          jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)               yqsol(:knon) = qsol(ni(:knon))
348            END DO            ELSE
349            CALL histend(nidbg)               yqsol = 0.
350            CALL histsync(nidbg)            END IF
351        ENDIF  
352                        DO k = 1, nsoilmx
353        DO k = 1, klev   ! epaisseur de couche               DO j = 1, knon
354        DO i = 1, klon                  i = ni(j)
355           delp(i,k) = paprs(i,k)-paprs(i,k+1)                  ytsoil(j, k) = ftsoil(i, k, nsrf)
356        ENDDO               END DO
357        ENDDO            END DO
358        DO i = 1, klon  ! vent de la premiere couche  
359           zx_alf1 = 1.0            DO k = 1, klev
360           zx_alf2 = 1.0 - zx_alf1               DO j = 1, knon
361           u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2                  i = ni(j)
362           v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2                  ypaprs(j, k) = paprs(i, k)
363        ENDDO                  ypplay(j, k) = pplay(i, k)
364  c                  ydelp(j, k) = delp(i, k)
365  c initialisation:                  yu(j, k) = u(i, k)
366  c                  yv(j, k) = v(i, k)
367        DO i = 1, klon                  yt(j, k) = t(i, k)
368           rugmer(i) = 0.0                  yq(j, k) = q(i, k)
369           cdragh(i) = 0.0               END DO
370           cdragm(i) = 0.0            END DO
371           dflux_t(i) = 0.0  
372           dflux_q(i) = 0.0            ! calculer Cdrag et les coefficients d'echange
373           zu1(i) = 0.0            CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &
374           zv1(i) = 0.0                 yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))
375        ENDDO            IF (iflag_pbl == 1) THEN
376        ypct = 0.0               CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
377        yts = 0.0               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
378        ysnow = 0.0               coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
379        yqsurf = 0.0            END IF
380        yalb = 0.0  
381        yalblw = 0.0            ! on met un seuil pour coefm et coefh
382        yrain_f = 0.0            IF (nsrf == is_oce) THEN
383        ysnow_f = 0.0               coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
384        yfder = 0.0               coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
385        ytaux = 0.0            END IF
386        ytauy = 0.0  
387        ysolsw = 0.0            IF (ok_kzmin) THEN
388        ysollw = 0.0               ! Calcul d'une diffusion minimale pour les conditions tres stables
389        ysollwdown = 0.0               CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
390        yrugos = 0.0                    coefm(:knon, 1), ycoefm0, ycoefh0)
391        yu1 = 0.0               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
392        yv1 = 0.0               coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
393        yrads = 0.0            END IF
394        ypaprs = 0.0  
395        ypplay = 0.0            IF (iflag_pbl >= 3) THEN
396        ydelp = 0.0               ! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et
397        yu = 0.0               ! Fr\'ed\'eric Hourdin
398        yv = 0.0               yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
399        yt = 0.0                    + ypplay(:knon, 1))) &
400        yq = 0.0                    * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
401        pctsrf_new = 0.0               DO k = 2, klev
402        y_flux_u = 0.0                  yzlay(1:knon, k) = yzlay(1:knon, k-1) &
403        y_flux_v = 0.0                       + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
404  C$$ PB                       / ypaprs(1:knon, k) &
405        y_dflux_t = 0.0                       * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
406        y_dflux_q = 0.0               END DO
407        ytsoil = 999999.               DO k = 1, klev
408        yrugoro = 0.                  yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
409  c -- LOOP                       / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
410        yu10mx = 0.0               END DO
411        yu10my = 0.0               yzlev(1:knon, 1) = 0.
412        ywindsp = 0.0               yzlev(:knon, klev+1) = 2. * yzlay(:knon, klev) &
413  c -- LOOP                    - yzlay(:knon, klev - 1)
414        DO nsrf = 1, nbsrf               DO k = 2, klev
415        DO i = 1, klon                  yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
416           d_ts(i,nsrf) = 0.0               END DO
417        ENDDO               DO k = 1, klev + 1
418        END DO                  DO j = 1, knon
419  C§§§ PB                     i = ni(j)
420        yfluxlat=0.                     yq2(j, k) = q2(i, k, nsrf)
421        flux_t = 0.                  END DO
422        flux_q = 0.               END DO
423        flux_u = 0.  
424        flux_v = 0.               CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
425        DO k = 1, klev               IF (prt_level > 9) PRINT *, 'USTAR = ', yustar
426        DO i = 1, klon  
427           d_t(i,k) = 0.0               ! iflag_pbl peut \^etre utilis\'e comme longueur de m\'elange
428           d_q(i,k) = 0.0  
429  c$$$         flux_t(i,k) = 0.0               IF (iflag_pbl >= 11) THEN
430  c$$$         flux_q(i,k) = 0.0                  CALL vdif_kcay(knon, dtime, rg, ypaprs, yzlev, yzlay, yu, yv, &
431           d_u(i,k) = 0.0                       yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, yustar, &
432           d_v(i,k) = 0.0                       iflag_pbl)
433  c$$$         flux_u(i,k) = 0.0               ELSE
434  c$$$         flux_v(i,k) = 0.0                  CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
435           zcoefh(i,k) = 0.0                       coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
436        ENDDO               END IF
437        ENDDO  
438  cAA      IF (itr.GE.1) THEN               coefm(:knon, 2:) = ykmm(:knon, 2:klev)
439  cAA      DO it = 1, itr               coefh(:knon, 2:) = ykmn(:knon, 2:klev)
440  cAA      DO k = 1, klev            END IF
441  cAA      DO i = 1, klon  
442  cAA         d_tr(i,k,it) = 0.0            ! calculer la diffusion des vitesses "u" et "v"
443  cAA      ENDDO            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &
444  cAA      ENDDO                 ypplay, ydelp, y_d_u, y_flux_u(:knon))
445  cAA      ENDDO            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
446  cAA      ENDIF                 ypplay, ydelp, y_d_v, y_flux_v(:knon))
447    
448  c            ! calculer la diffusion de "q" et de "h"
449  c Boucler sur toutes les sous-fractions du sol:            CALL clqh(dtime, jour, firstcal, rlat, nsrf, ni(:knon), ytsoil, &
450  c                 yqsol, rmu0, yrugos, yrugoro, yu1, yv1, coefh(:knon, :), yt, &
451  C Initialisation des "pourcentages potentiels". On considere ici qu'on                 yq, yts(:knon), ypaprs, ypplay, ydelp, yrads, yalb(:knon), &
452  C peut avoir potentiellementdela glace sur tout le domaine oceanique                 ysnow, yqsurf, yrain_f, ysnow_f, yfder, yfluxlat, &
453  C (a affiner)                 pctsrf_new_sic, yagesno(:knon), y_d_t, y_d_q, y_d_ts(:knon), &
454                   yz0_new, y_flux_t(:knon), y_flux_q(:knon), y_dflux_t, &
455        pctsrf_pot = pctsrf                 y_dflux_q, y_fqcalving, y_ffonte, y_run_off_lic_0)
456        pctsrf_pot(:,is_oce) = 1. - zmasq(:)  
457        pctsrf_pot(:,is_sic) = 1. - zmasq(:)            ! calculer la longueur de rugosite sur ocean
458              yrugm = 0.
459        DO 99999 nsrf = 1, nbsrf            IF (nsrf == is_oce) THEN
460                 DO j = 1, knon
461  c chercher les indices:                  yrugm(j) = 0.018*coefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
462        DO j = 1, klon                       0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))
463           ni(j) = 0                  yrugm(j) = max(1.5E-05, yrugm(j))
464        ENDDO               END DO
465        knon = 0            END IF
466        DO i = 1, klon            DO j = 1, knon
467                 y_dflux_t(j) = y_dflux_t(j)*ypct(j)
468  C pour determiner le domaine a traiter on utilise les surfaces "potentielles"               y_dflux_q(j) = y_dflux_q(j)*ypct(j)
469  C                 yu1(j) = yu1(j)*ypct(j)
470        IF (pctsrf_pot(i,nsrf).GT.epsfra) THEN               yv1(j) = yv1(j)*ypct(j)
471           knon = knon + 1            END DO
472           ni(knon) = i  
473        ENDIF            DO k = 1, klev
474        ENDDO               DO j = 1, knon
475  c                  i = ni(j)
476        if (check) THEN                  coefh(j, k) = coefh(j, k)*ypct(j)
477            write(*,*)'CLMAIN, nsrf, knon =',nsrf, knon                  coefm(j, k) = coefm(j, k)*ypct(j)
478  CC        call flush(6)                  y_d_t(j, k) = y_d_t(j, k)*ypct(j)
479        endif                  y_d_q(j, k) = y_d_q(j, k)*ypct(j)
480  c                  y_d_u(j, k) = y_d_u(j, k)*ypct(j)
481  c variables pour avoir une sortie IOIPSL des INDEX                  y_d_v(j, k) = y_d_v(j, k)*ypct(j)
482  c               END DO
483        IF (debugindex) THEN            END DO
484            tabindx(:)=0.  
485  c          tabindx(1:knon)=(/FLOAT(i),i=1:knon/)            DO j = 1, knon
486            DO i=1,knon               i = ni(j)
487              tabindx(1:knon)=FLOAT(i)               flux_t(i, nsrf) = y_flux_t(j)
488            END DO               flux_q(i, nsrf) = y_flux_q(j)
489            debugtab(:,:)=0.               flux_u(i, nsrf) = y_flux_u(j)
490            ndexbg(:)=0               flux_v(i, nsrf) = y_flux_v(j)
491            CALL gath2cpl(tabindx,debugtab,klon,knon,iim,jjm,ni)            END DO
492            CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1)  
493       $        ,ndexbg)            evap(:, nsrf) = -flux_q(:, nsrf)
494        ENDIF  
495        IF (knon.EQ.0) GOTO 99999            falbe(:, nsrf) = 0.
496        DO j = 1, knon            snow(:, nsrf) = 0.
497        i = ni(j)            qsurf(:, nsrf) = 0.
498          ypct(j) = pctsrf(i,nsrf)            rugos(:, nsrf) = 0.
499          yts(j) = ts(i,nsrf)            fluxlat(:, nsrf) = 0.
500  cIM "slab" ocean            DO j = 1, knon
501  c        PRINT *, 'tslab = ', i, tslab(i)               i = ni(j)
502          ytslab(i) = tslab(i)               d_ts(i, nsrf) = y_d_ts(j)
503  c               falbe(i, nsrf) = yalb(j)
504          ysnow(j) = snow(i,nsrf)               snow(i, nsrf) = ysnow(j)
505          yqsurf(j) = qsurf(i,nsrf)               qsurf(i, nsrf) = yqsurf(j)
506          yalb(j) = albe(i,nsrf)               rugos(i, nsrf) = yz0_new(j)
507          yalblw(j) = alblw(i,nsrf)               fluxlat(i, nsrf) = yfluxlat(j)
508          yrain_f(j) = rain_f(i)               IF (nsrf == is_oce) THEN
509          ysnow_f(j) = snow_f(i)                  rugmer(i) = yrugm(j)
510          yagesno(j) = agesno(i,nsrf)                  rugos(i, nsrf) = yrugm(j)
511          yfder(j) = fder(i)               END IF
512          ytaux(j) = flux_u(i,1,nsrf)               agesno(i, nsrf) = yagesno(j)
513          ytauy(j) = flux_v(i,1,nsrf)               fqcalving(i, nsrf) = y_fqcalving(j)
514          ysolsw(j) = solsw(i,nsrf)               ffonte(i, nsrf) = y_ffonte(j)
515          ysollw(j) = sollw(i,nsrf)               cdragh(i) = cdragh(i) + coefh(j, 1)
516          ysollwdown(j) = sollwdown(i)               cdragm(i) = cdragm(i) + coefm(j, 1)
517          yrugos(j) = rugos(i,nsrf)               dflux_t(i) = dflux_t(i) + y_dflux_t(j)
518          yrugoro(j) = rugoro(i)               dflux_q(i) = dflux_q(i) + y_dflux_q(j)
519          yu1(j) = u1lay(i)               zu1(i) = zu1(i) + yu1(j)
520          yv1(j) = v1lay(i)               zv1(i) = zv1(i) + yv1(j)
521          yrads(j) =  ysolsw(j)+ ysollw(j)            END DO
522          ypaprs(j,klev+1) = paprs(i,klev+1)            IF (nsrf == is_ter) THEN
523          y_run_off_lic_0(j) = run_off_lic_0(i)               qsol(ni(:knon)) = yqsol(:knon)
524  c -- LOOP            else IF (nsrf == is_lic) THEN
525         yu10mx(j) = u10m(i,nsrf)               DO j = 1, knon
526         yu10my(j) = v10m(i,nsrf)                  i = ni(j)
527         ywindsp(j) = SQRT(yu10mx(j)*yu10mx(j) + yu10my(j)*yu10my(j) )                  run_off_lic_0(i) = y_run_off_lic_0(j)
528  c -- LOOP               END DO
529        END DO            END IF
530  C  
531  C     IF bucket model for continent, copy soil water content            ftsoil(:, :, nsrf) = 0.
532        IF ( nsrf .eq. is_ter .and. .not. ok_veget ) THEN            DO k = 1, nsoilmx
533                 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
540                 i = ni(j)
541                 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
549    
550              ! diagnostic t, q a 2m et u, v a 10m
551    
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            DO j = 1, knon
576              i = ni(j)               i = ni(j)
577              yqsol(j) = qsol(i)               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(ypaprs, ypplay, yt2m, yq2m, yustar, y_flux_t(:knon), &
587                   y_flux_q(:knon), yu, yv, yt, yq, ypblh(:knon), ycapcl, &
588                   yoliqcl, 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            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  
603    
604        RETURN            DO j = 1, knon
605        END               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.10  
changed lines
  Added in v.207

  ViewVC Help
Powered by ViewVC 1.1.21