/[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 191 by guez, Mon May 9 19:56:28 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, pctsrf_new, t, q, u, v, jour, rmu0, ts, &
8       .                  fluxlat,         cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, pplay, &
9       .                  rain_f, snow_f, solsw, sollw, sollwdown, fder,         snow, qsurf, evap, falbe, fluxlat, rain_fall, snow_f, solsw, sollw, &
10       .                  rlon, rlat, cufi, cvfi, rugos,         fder, rlat, rugos, firstcal, agesno, rugoro, d_t, d_q, d_u, d_v, d_ts, &
11       .                  debut, lafin, agesno,rugoro,         flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, dflux_t, dflux_q, &
12       .                  d_t,d_q,d_u,d_v,d_ts,         ycoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh, capcl, oliqcl, cteicl, &
13       .                  flux_t,flux_q,flux_u,flux_v,cdragh,cdragm,         pblt, therm, 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
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 stdlevvar_m, only: stdlevvar
40  cAA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir      USE suphec_m, ONLY: rd, rg, rkappa
41  cAA si les informations des subsurfaces doivent etre prises en compte      use ustarhb_m, only: ustarhb
42  cAA il faudra sortir ces memes champs en leur ajoutant une dimension,      use vdif_kcay_m, only: vdif_kcay
43  cAA c'est a dire nbsrf (nbre de subsurface).      use yamada4_m, only: yamada4
44        USE ioipsl  
45        USE interface_surf      REAL, INTENT(IN):: dtime ! interval du temps (secondes)
46        use dimens_m      REAL, INTENT(inout):: pctsrf(klon, nbsrf)
47        use indicesol  
48        use dimphy      ! la nouvelle repartition des surfaces sortie de l'interface
49        use dimsoil      REAL, INTENT(out):: pctsrf_new(klon, nbsrf)
50        use temps  
51        use iniprint      REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
52        use YOMCST      REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg/kg)
53        use yoethf      REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
54        use fcttre      INTEGER, INTENT(IN):: jour ! jour de l'annee en cours
55        use conf_phys_m      REAL, intent(in):: rmu0(klon) ! cosinus de l'angle solaire zenithal    
56        use gath_cpl, only: gath2cpl      REAL, INTENT(IN):: ts(klon, nbsrf) ! temperature du sol (en Kelvin)
57        IMPLICIT none      REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
58  c======================================================================      REAL, INTENT(IN):: ksta, ksta_ter
59  c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818      LOGICAL, INTENT(IN):: ok_kzmin
60  c Objet: interface de "couche limite" (diffusion verticale)  
61  c Arguments:      REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
62  c dtime----input-R- interval du temps (secondes)      ! soil temperature of surface fraction
63  c itap-----input-I- numero du pas de temps  
64  c date0----input-R- jour initial      REAL, INTENT(inout):: qsol(klon)
65  c t--------input-R- temperature (K)      ! column-density of water in soil, in kg m-2
66  c q--------input-R- vapeur d'eau (kg/kg)  
67  c u--------input-R- vitesse u      REAL, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)
68  c v--------input-R- vitesse v      REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
69  c ts-------input-R- temperature du sol (en Kelvin)      REAL, INTENT(inout):: snow(klon, nbsrf)
70  c paprs----input-R- pression a intercouche (Pa)      REAL qsurf(klon, nbsrf)
71  c pplay----input-R- pression au milieu de couche (Pa)      REAL evap(klon, nbsrf)
72  c radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2      REAL, intent(inout):: falbe(klon, nbsrf)
73  c rlat-----input-R- latitude en degree  
74  c rugos----input-R- longeur de rugosite (en m)      REAL fluxlat(klon, nbsrf)
75  c cufi-----input-R- resolution des mailles en x (m)  
76  c cvfi-----input-R- resolution des mailles en y (m)      REAL, intent(in):: rain_fall(klon)
77  c      ! liquid water mass flux (kg/m2/s), positive down
78  c d_t------output-R- le changement pour "t"  
79  c d_q------output-R- le changement pour "q"      REAL, intent(in):: snow_f(klon)
80  c d_u------output-R- le changement pour "u"      ! solid water mass flux (kg/m2/s), positive down
81  c d_v------output-R- le changement pour "v"  
82  c d_ts-----output-R- le changement pour "ts"      REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)
83  c flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)      REAL, intent(in):: fder(klon)
84  c                    (orientation positive vers le bas)      REAL, INTENT(IN):: rlat(klon) ! latitude en degr\'es
85  c flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)  
86  c flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal      REAL, intent(inout):: rugos(klon, nbsrf) ! longueur de rugosit\'e (en m)
87  c flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal  
88  c dflux_t derive du flux sensible      LOGICAL, INTENT(IN):: firstcal
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      REAL y_fqcalving(klon), y_ffonte(klon)
150        logical, intent(in):: lafin      real y_run_off_lic_0(klon)
151        logical ok_veget  
152        character*6 ocean      REAL rugmer(klon)
153        integer npas, nexca  
154  c      REAL ytsoil(klon, nsoilmx)
155        REAL pctsrf(klon,nbsrf)  
156        REAL ts(klon,nbsrf)      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
157        REAL d_ts(klon,nbsrf)      REAL yalb(klon)
158        REAL snow(klon,nbsrf)      REAL yu1(klon), yv1(klon)
159        REAL qsurf(klon,nbsrf)      ! on rajoute en output yu1 et yv1 qui sont les vents dans
160        REAL evap(klon,nbsrf)      ! la premiere couche
161        REAL albe(klon,nbsrf)      REAL ysnow(klon), yqsurf(klon), yagesno(klon)
162        REAL alblw(klon,nbsrf)  
163  c$$$ PB      real yqsol(klon)
164        REAL fluxlat(klon,nbsrf)      ! column-density of water in soil, in kg m-2
165  C  
166        real rain_f(klon), snow_f(klon)      REAL yrain_f(klon)
167        REAL fder(klon)      ! liquid water mass flux (kg/m2/s), positive down
168  cIM cf. JLD   REAL sollw(klon), solsw(klon), sollwdown(klon)  
169        REAL sollw(klon,nbsrf), solsw(klon,nbsrf), sollwdown(klon)      REAL ysnow_f(klon)
170        REAL rugos(klon,nbsrf)      ! solid water mass flux (kg/m2/s), positive down
171  C la nouvelle repartition des surfaces sortie de l'interface  
172        REAL pctsrf_new(klon,nbsrf)      REAL yfder(klon)
173  cAA      REAL yrugm(klon), yrads(klon), yrugoro(klon)
174        REAL zcoefh(klon,klev)  
175        REAL zu1(klon)      REAL yfluxlat(klon)
176        REAL zv1(klon)  
177  cAA      REAL y_d_ts(klon)
178  c$$$ PB ajout pour soil      REAL y_d_t(klon, klev), y_d_q(klon, klev)
179        LOGICAL soil_model      REAL y_d_u(klon, klev), y_d_v(klon, klev)
180  cIM ajout seuils cdrm, cdrh      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)
181        REAL cdmmax, cdhmax      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)
182  cIM: 261103      REAL y_dflux_t(klon), y_dflux_q(klon)
183        REAL ksta, ksta_ter      REAL coefh(klon, klev), coefm(klon, klev)
184        LOGICAL ok_kzmin      REAL yu(klon, klev), yv(klon, klev)
185  cIM: 261103      REAL yt(klon, klev), yq(klon, klev)
186        REAL ftsoil(klon,nsoilmx,nbsrf)      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)
187        REAL ytsoil(klon,nsoilmx)  
188        REAL qsol(klon)      REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
189  c======================================================================  
190        EXTERNAL clqh, clvent, coefkz, calbeta, cltrac      REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
191  c======================================================================      REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
192        REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      REAL ykmq(klon, klev+1)
193        REAL yalb(klon)      REAL yq2(klon, klev+1)
194        REAL yalblw(klon)      REAL q2diag(klon, klev+1)
195        REAL yu1(klon), yv1(klon)  
196        real ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)      REAL u1lay(klon), v1lay(klon)
197        real yrain_f(klon), ysnow_f(klon)      REAL delp(klon, klev)
198        real ysollw(klon), ysolsw(klon), ysollwdown(klon)      INTEGER i, k, nsrf
199        real yfder(klon), ytaux(klon), ytauy(klon)  
200        REAL yrugm(klon), yrads(klon),yrugoro(klon)      INTEGER ni(klon), knon, j
201  c$$$ PB  
202        REAL yfluxlat(klon)      REAL pctsrf_pot(klon, nbsrf)
203  C      ! "pourcentage potentiel" pour tenir compte des \'eventuelles
204        REAL y_d_ts(klon)      ! apparitions ou disparitions de la glace de mer
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 zx_alf1, zx_alf2 ! valeur ambiante par extrapolation
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 yt2m(klon), yq2m(klon), yu10m(klon)
209        REAL y_dflux_t(klon), y_dflux_q(klon)      REAL yustar(klon)
210        REAL ycoefh(klon,klev), ycoefm(klon,klev)  
211        REAL yu(klon,klev), yv(klon,klev)      REAL yt10m(klon), yq10m(klon)
212        REAL yt(klon,klev), yq(klon,klev)      REAL ypblh(klon)
213        REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)      REAL ylcl(klon)
214  c      REAL ycapcl(klon)
215        LOGICAL ok_nonloc      REAL yoliqcl(klon)
216        PARAMETER (ok_nonloc=.FALSE.)      REAL ycteicl(klon)
217        REAL ycoefm0(klon,klev), ycoefh0(klon,klev)      REAL ypblt(klon)
218        REAL ytherm(klon)
219  cIM 081204 hcl_Anne ? BEG      REAL ytrmb1(klon)
220        real yzlay(klon,klev),yzlev(klon,klev+1),yteta(klon,klev)      REAL ytrmb2(klon)
221        real ykmm(klon,klev+1),ykmn(klon,klev+1)      REAL ytrmb3(klon)
222        real ykmq(klon,klev+1)      REAL uzon(klon), vmer(klon)
223        real yq2(klon,klev+1),q2(klon,klev+1,nbsrf)      REAL tair1(klon), qair1(klon), tairsol(klon)
224        real q2diag(klon,klev+1)      REAL psfce(klon), patm(klon)
225  cIM 081204   real yustar(klon),y_cd_m(klon),y_cd_h(klon)  
226  cIM 081204 hcl_Anne ? END      REAL qairsol(klon), zgeo1(klon)
227  c      REAL rugo1(klon)
228        REAL u1lay(klon), v1lay(klon)  
229        REAL delp(klon,klev)      ! utiliser un jeu de fonctions simples              
230        INTEGER i, k, nsrf      LOGICAL zxli
231  cAA   INTEGER it      PARAMETER (zxli=.FALSE.)
232        INTEGER ni(klon), knon, j  
233  c Introduction d'une variable "pourcentage potentiel" pour tenir compte      !------------------------------------------------------------
234  c des eventuelles apparitions et/ou disparitions de la glace de mer  
235        REAL pctsrf_pot(klon,nbsrf)      ytherm = 0.
236          
237  c======================================================================      DO k = 1, klev ! epaisseur de couche
238        REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.         DO i = 1, klon
239  c======================================================================            delp(i, k) = paprs(i, k) - paprs(i, k+1)
240  c         END DO
241  c maf pour sorties IOISPL en cas de debugagage      END DO
242  c      DO i = 1, klon ! vent de la premiere couche
243        CHARACTER*80 cldebug         zx_alf1 = 1.0
244        SAVE cldebug         zx_alf2 = 1.0 - zx_alf1
245        CHARACTER*8 cl_surf(nbsrf)         u1lay(i) = u(i, 1)*zx_alf1 + u(i, 2)*zx_alf2
246        SAVE cl_surf         v1lay(i) = v(i, 1)*zx_alf1 + v(i, 2)*zx_alf2
247        INTEGER nhoridbg, nidbg      END DO
248        SAVE nhoridbg, nidbg  
249        INTEGER ndexbg(iim*(jjm+1))      ! Initialization:
250        REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1), zjulian      rugmer = 0.
251        REAL tabindx(klon)      cdragh = 0.
252        REAL debugtab(iim,jjm+1)      cdragm = 0.
253        LOGICAL first_appel      dflux_t = 0.
254        SAVE first_appel      dflux_q = 0.
255        DATA first_appel/.true./      zu1 = 0.
256        LOGICAL debugindex      zv1 = 0.
257        SAVE debugindex      ypct = 0.
258        DATA debugindex/.false./      yts = 0.
259        integer idayref      ysnow = 0.
260        REAL t2m(klon,nbsrf), q2m(klon,nbsrf)      yqsurf = 0.
261        REAL u10m(klon,nbsrf), v10m(klon,nbsrf)      yrain_f = 0.
262  c      ysnow_f = 0.
263        REAL yt2m(klon), yq2m(klon), yu10m(klon)      yfder = 0.
264        REAL yustar(klon)      yrugos = 0.
265  c -- LOOP      yu1 = 0.
266         REAL yu10mx(klon)      yv1 = 0.
267         REAL yu10my(klon)      yrads = 0.
268         REAL ywindsp(klon)      ypaprs = 0.
269  c -- LOOP      ypplay = 0.
270  c      ydelp = 0.
271        REAL yt10m(klon), yq10m(klon)      yu = 0.
272  cIM cf. AM : pbl, hbtm2 (Comme les autres diagnostics on cumule ds physic ce qui      yv = 0.
273  c   permet de sortir les grdeurs par sous surface)      yt = 0.
274        REAL pblh(klon,nbsrf)      yq = 0.
275        REAL plcl(klon,nbsrf)      pctsrf_new = 0.
276        REAL capCL(klon,nbsrf)      y_flux_u = 0.
277        REAL oliqCL(klon,nbsrf)      y_flux_v = 0.
278        REAL cteiCL(klon,nbsrf)      y_dflux_t = 0.
279        REAL pblT(klon,nbsrf)      y_dflux_q = 0.
280        REAL therm(klon,nbsrf)      ytsoil = 999999.
281        REAL trmb1(klon,nbsrf)      yrugoro = 0.
282        REAL trmb2(klon,nbsrf)      d_ts = 0.
283        REAL trmb3(klon,nbsrf)      yfluxlat = 0.
284        REAL ypblh(klon)      flux_t = 0.
285        REAL ylcl(klon)      flux_q = 0.
286        REAL ycapCL(klon)      flux_u = 0.
287        REAL yoliqCL(klon)      flux_v = 0.
288        REAL ycteiCL(klon)      d_t = 0.
289        REAL ypblT(klon)      d_q = 0.
290        REAL ytherm(klon)      d_u = 0.
291        REAL ytrmb1(klon)      d_v = 0.
292        REAL ytrmb2(klon)      ycoefh = 0.
293        REAL ytrmb3(klon)  
294        REAL y_cd_h(klon), y_cd_m(klon)      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
295  c     REAL ygamt(klon,2:klev) ! contre-gradient pour temperature      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
296  c     REAL ygamq(klon,2:klev) ! contre-gradient pour humidite      ! (\`a affiner)
297        REAL uzon(klon), vmer(klon)  
298        REAL tair1(klon), qair1(klon), tairsol(klon)      pctsrf_pot = pctsrf
299        REAL psfce(klon), patm(klon)      pctsrf_pot(:, is_oce) = 1. - zmasq
300  c      pctsrf_pot(:, is_sic) = 1. - zmasq
301        REAL qairsol(klon), zgeo1(klon)  
302        REAL rugo1(klon)      ! Boucler sur toutes les sous-fractions du sol:
303  c  
304        LOGICAL zxli ! utiliser un jeu de fonctions simples      loop_surface: DO nsrf = 1, nbsrf
305        PARAMETER (zxli=.FALSE.)         ! Chercher les indices :
306  c         ni = 0
307        REAL zt, zqs, zdelta, zcor         knon = 0
308        REAL t_coup         DO i = 1, klon
309        PARAMETER(t_coup=273.15)            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
310  C            ! "potentielles"
311        character (len = 20) :: modname = 'clmain'            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
312        LOGICAL check               knon = knon + 1
313        PARAMETER (check=.false.)               ni(knon) = i
314              END IF
315           END DO
316  c initialisation Anne  
317        ytherm(:) = 0.         if_knon: IF (knon /= 0) then
318  C            DO j = 1, knon
319        if (check) THEN               i = ni(j)
320            write(*,*) modname,'  klon=',klon               ypct(j) = pctsrf(i, nsrf)
321  CC        call flush(6)               yts(j) = ts(i, nsrf)
322        endif               ysnow(j) = snow(i, nsrf)
323        IF (debugindex .and. first_appel) THEN               yqsurf(j) = qsurf(i, nsrf)
324            first_appel=.false.               yalb(j) = falbe(i, nsrf)
325  !               yrain_f(j) = rain_fall(i)
326  ! initialisation sorties netcdf               ysnow_f(j) = snow_f(i)
327  !               yagesno(j) = agesno(i, nsrf)
328            idayref = day_ini               yfder(j) = fder(i)
329            CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)               yrugos(j) = rugos(i, nsrf)
330            CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)               yrugoro(j) = rugoro(i)
331            DO i = 1, iim               yu1(j) = u1lay(i)
332              zx_lon(i,1) = rlon(i+1)               yv1(j) = v1lay(i)
333              zx_lon(i,jjm+1) = rlon(i+1)               yrads(j) = solsw(i, nsrf) + sollw(i, nsrf)
334            ENDDO               ypaprs(j, klev+1) = paprs(i, klev+1)
335            CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)               y_run_off_lic_0(j) = run_off_lic_0(i)
336            cldebug='sous_index'            END DO
337            CALL histbeg_totreg(cldebug, iim,zx_lon(:,1),jjm+1,  
338       $        zx_lat(1,:),1,iim,1,jjm            ! For continent, copy soil water content
339       $        +1, itau_phy,zjulian,dtime,nhoridbg,nidbg)            IF (nsrf == is_ter) THEN
340  ! no vertical axis               yqsol(:knon) = qsol(ni(:knon))
341            cl_surf(1)='ter'            ELSE
342            cl_surf(2)='lic'               yqsol = 0.
343            cl_surf(3)='oce'            END IF
344            cl_surf(4)='sic'  
345            DO nsrf=1,nbsrf            DO k = 1, nsoilmx
346              CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim,               DO j = 1, knon
347       $          jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)                  i = ni(j)
348            END DO                  ytsoil(j, k) = ftsoil(i, k, nsrf)
349            CALL histend(nidbg)               END DO
350            CALL histsync(nidbg)            END DO
351        ENDIF  
352                        DO k = 1, klev
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)                  ypaprs(j, k) = paprs(i, k)
356        ENDDO                  ypplay(j, k) = pplay(i, k)
357        ENDDO                  ydelp(j, k) = delp(i, k)
358        DO i = 1, klon  ! vent de la premiere couche                  yu(j, k) = u(i, k)
359  ccc         zx_alf1 = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))                  yv(j, k) = v(i, k)
360           zx_alf1 = 1.0                  yt(j, k) = t(i, k)
361           zx_alf2 = 1.0 - zx_alf1                  yq(j, k) = q(i, k)
362           u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2               END DO
363           v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2            END DO
364        ENDDO  
365  c            ! calculer Cdrag et les coefficients d'echange
366  c initialisation:            CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &
367  c                 yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))
368        DO i = 1, klon            IF (iflag_pbl == 1) THEN
369           rugmer(i) = 0.0               CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
370           cdragh(i) = 0.0               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
371           cdragm(i) = 0.0               coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
372           dflux_t(i) = 0.0            END IF
373           dflux_q(i) = 0.0  
374           zu1(i) = 0.0            ! on met un seuil pour coefm et coefh
375           zv1(i) = 0.0            IF (nsrf == is_oce) THEN
376        ENDDO               coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
377        ypct = 0.0               coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
378        yts = 0.0            END IF
379        ysnow = 0.0  
380        yqsurf = 0.0            IF (ok_kzmin) THEN
381        yalb = 0.0               ! Calcul d'une diffusion minimale pour les conditions tres stables
382        yalblw = 0.0               CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
383        yrain_f = 0.0                    coefm(:knon, 1), ycoefm0, ycoefh0)
384        ysnow_f = 0.0               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
385        yfder = 0.0               coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
386        ytaux = 0.0            END IF
387        ytauy = 0.0  
388        ysolsw = 0.0            IF (iflag_pbl >= 3) THEN
389        ysollw = 0.0               ! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et
390        ysollwdown = 0.0               ! Fr\'ed\'eric Hourdin
391        yrugos = 0.0               yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
392        yu1 = 0.0                    + ypplay(:knon, 1))) &
393        yv1 = 0.0                    * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
394        yrads = 0.0               DO k = 2, klev
395        ypaprs = 0.0                  yzlay(1:knon, k) = yzlay(1:knon, k-1) &
396        ypplay = 0.0                       + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
397        ydelp = 0.0                       / ypaprs(1:knon, k) &
398        yu = 0.0                       * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
399        yv = 0.0               END DO
400        yt = 0.0               DO k = 1, klev
401        yq = 0.0                  yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
402        pctsrf_new = 0.0                       / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
403        y_flux_u = 0.0               END DO
404        y_flux_v = 0.0               yzlev(1:knon, 1) = 0.
405  C$$ PB               yzlev(:knon, klev+1) = 2. * yzlay(:knon, klev) &
406        y_dflux_t = 0.0                    - yzlay(:knon, klev - 1)
407        y_dflux_q = 0.0               DO k = 2, klev
408        ytsoil = 999999.                  yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
409        yrugoro = 0.               END DO
410  c -- LOOP               DO k = 1, klev + 1
411        yu10mx = 0.0                  DO j = 1, knon
412        yu10my = 0.0                     i = ni(j)
413        ywindsp = 0.0                     yq2(j, k) = q2(i, k, nsrf)
414  c -- LOOP                  END DO
415        DO nsrf = 1, nbsrf               END DO
416        DO i = 1, klon  
417           d_ts(i,nsrf) = 0.0               CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
418        ENDDO               IF (prt_level > 9) PRINT *, 'USTAR = ', yustar
419        END DO  
420  C§§§ PB               ! iflag_pbl peut \^etre utilis\'e comme longueur de m\'elange
421        yfluxlat=0.  
422        flux_t = 0.               IF (iflag_pbl >= 11) THEN
423        flux_q = 0.                  CALL vdif_kcay(knon, dtime, rg, ypaprs, yzlev, yzlay, yu, yv, &
424        flux_u = 0.                       yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, yustar, &
425        flux_v = 0.                       iflag_pbl)
426        DO k = 1, klev               ELSE
427        DO i = 1, klon                  CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
428           d_t(i,k) = 0.0                       coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
429           d_q(i,k) = 0.0               END IF
430  c$$$         flux_t(i,k) = 0.0  
431  c$$$         flux_q(i,k) = 0.0               coefm(:knon, 2:) = ykmm(:knon, 2:klev)
432           d_u(i,k) = 0.0               coefh(:knon, 2:) = ykmn(:knon, 2:klev)
433           d_v(i,k) = 0.0            END IF
434  c$$$         flux_u(i,k) = 0.0  
435  c$$$         flux_v(i,k) = 0.0            ! calculer la diffusion des vitesses "u" et "v"
436           zcoefh(i,k) = 0.0            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &
437        ENDDO                 ypplay, ydelp, y_d_u, y_flux_u)
438        ENDDO            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
439  cAA      IF (itr.GE.1) THEN                 ypplay, ydelp, y_d_v, y_flux_v)
440  cAA      DO it = 1, itr  
441  cAA      DO k = 1, klev            ! calculer la diffusion de "q" et de "h"
442  cAA      DO i = 1, klon            CALL clqh(dtime, jour, firstcal, rlat, knon, nsrf, ni(:knon), &
443  cAA         d_tr(i,k,it) = 0.0                 pctsrf, ytsoil, yqsol, rmu0, yrugos, yrugoro, yu1, yv1, &
444  cAA      ENDDO                 coefh(:knon, :), yt, yq, yts, ypaprs, ypplay, ydelp, yrads, &
445  cAA      ENDDO                 yalb(:knon), ysnow, yqsurf, yrain_f, ysnow_f, yfder, yfluxlat, &
446  cAA      ENDDO                 pctsrf_new, yagesno(:knon), y_d_t, y_d_q, y_d_ts(:knon), &
447  cAA      ENDIF                 yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q, &
448                   y_fqcalving, y_ffonte, y_run_off_lic_0)
449  c  
450  c Boucler sur toutes les sous-fractions du sol:            ! calculer la longueur de rugosite sur ocean
451  c            yrugm = 0.
452  C Initialisation des "pourcentages potentiels". On considere ici qu'on            IF (nsrf == is_oce) THEN
453  C peut avoir potentiellementdela glace sur tout le domaine oceanique               DO j = 1, knon
454  C (a affiner)                  yrugm(j) = 0.018*coefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
455                         0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))
456        pctsrf_pot = pctsrf                  yrugm(j) = max(1.5E-05, yrugm(j))
457        pctsrf_pot(:,is_oce) = 1. - zmasq(:)               END DO
458        pctsrf_pot(:,is_sic) = 1. - zmasq(:)            END IF
459              DO j = 1, knon
460        DO 99999 nsrf = 1, nbsrf               y_dflux_t(j) = y_dflux_t(j)*ypct(j)
461                 y_dflux_q(j) = y_dflux_q(j)*ypct(j)
462  c chercher les indices:               yu1(j) = yu1(j)*ypct(j)
463        DO j = 1, klon               yv1(j) = yv1(j)*ypct(j)
464           ni(j) = 0            END DO
465        ENDDO  
466        knon = 0            DO k = 1, klev
467        DO i = 1, klon               DO j = 1, knon
468                    i = ni(j)
469  C pour determiner le domaine a traiter on utilise les surfaces "potentielles"                  coefh(j, k) = coefh(j, k)*ypct(j)
470  C                    coefm(j, k) = coefm(j, k)*ypct(j)
471        IF (pctsrf_pot(i,nsrf).GT.epsfra) THEN                  y_d_t(j, k) = y_d_t(j, k)*ypct(j)
472           knon = knon + 1                  y_d_q(j, k) = y_d_q(j, k)*ypct(j)
473           ni(knon) = i                  flux_t(i, k, nsrf) = y_flux_t(j, k)
474        ENDIF                  flux_q(i, k, nsrf) = y_flux_q(j, k)
475        ENDDO                  flux_u(i, k, nsrf) = y_flux_u(j, k)
476  c                  flux_v(i, k, nsrf) = y_flux_v(j, k)
477        if (check) THEN                  y_d_u(j, k) = y_d_u(j, k)*ypct(j)
478            write(*,*)'CLMAIN, nsrf, knon =',nsrf, knon                  y_d_v(j, k) = y_d_v(j, k)*ypct(j)
479  CC        call flush(6)               END DO
480        endif            END DO
481  c  
482  c variables pour avoir une sortie IOIPSL des INDEX            evap(:, nsrf) = -flux_q(:, 1, nsrf)
483  c  
484        IF (debugindex) THEN            falbe(:, nsrf) = 0.
485            tabindx(:)=0.            snow(:, nsrf) = 0.
486  c          tabindx(1:knon)=(/FLOAT(i),i=1:knon/)            qsurf(:, nsrf) = 0.
487            DO i=1,knon            rugos(:, nsrf) = 0.
488              tabindx(1:knon)=FLOAT(i)            fluxlat(:, nsrf) = 0.
489            END DO            DO j = 1, knon
490            debugtab(:,:)=0.               i = ni(j)
491            ndexbg(:)=0               d_ts(i, nsrf) = y_d_ts(j)
492            CALL gath2cpl(tabindx,debugtab,klon,knon,iim,jjm,ni)               falbe(i, nsrf) = yalb(j)
493            CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1)               snow(i, nsrf) = ysnow(j)
494       $        ,ndexbg)               qsurf(i, nsrf) = yqsurf(j)
495        ENDIF               rugos(i, nsrf) = yz0_new(j)
496        IF (knon.EQ.0) GOTO 99999               fluxlat(i, nsrf) = yfluxlat(j)
497        DO j = 1, knon               IF (nsrf == is_oce) THEN
498        i = ni(j)                  rugmer(i) = yrugm(j)
499          ypct(j) = pctsrf(i,nsrf)                  rugos(i, nsrf) = yrugm(j)
500          yts(j) = ts(i,nsrf)               END IF
501  cIM "slab" ocean               agesno(i, nsrf) = yagesno(j)
502  c        PRINT *, 'tslab = ', i, tslab(i)               fqcalving(i, nsrf) = y_fqcalving(j)
503          ytslab(i) = tslab(i)               ffonte(i, nsrf) = y_ffonte(j)
504  c               cdragh(i) = cdragh(i) + coefh(j, 1)
505          ysnow(j) = snow(i,nsrf)               cdragm(i) = cdragm(i) + coefm(j, 1)
506          yqsurf(j) = qsurf(i,nsrf)               dflux_t(i) = dflux_t(i) + y_dflux_t(j)
507          yalb(j) = albe(i,nsrf)               dflux_q(i) = dflux_q(i) + y_dflux_q(j)
508          yalblw(j) = alblw(i,nsrf)               zu1(i) = zu1(i) + yu1(j)
509          yrain_f(j) = rain_f(i)               zv1(i) = zv1(i) + yv1(j)
510          ysnow_f(j) = snow_f(i)            END DO
511          yagesno(j) = agesno(i,nsrf)            IF (nsrf == is_ter) THEN
512          yfder(j) = fder(i)               qsol(ni(:knon)) = yqsol(:knon)
513          ytaux(j) = flux_u(i,1,nsrf)            else IF (nsrf == is_lic) THEN
514          ytauy(j) = flux_v(i,1,nsrf)               DO j = 1, knon
515          ysolsw(j) = solsw(i,nsrf)                  i = ni(j)
516          ysollw(j) = sollw(i,nsrf)                  run_off_lic_0(i) = y_run_off_lic_0(j)
517          ysollwdown(j) = sollwdown(i)               END DO
518          yrugos(j) = rugos(i,nsrf)            END IF
519          yrugoro(j) = rugoro(i)  
520          yu1(j) = u1lay(i)            ftsoil(:, :, nsrf) = 0.
521          yv1(j) = v1lay(i)            DO k = 1, nsoilmx
522          yrads(j) =  ysolsw(j)+ ysollw(j)               DO j = 1, knon
523          ypaprs(j,klev+1) = paprs(i,klev+1)                  i = ni(j)
524          y_run_off_lic_0(j) = run_off_lic_0(i)                  ftsoil(i, k, nsrf) = ytsoil(j, k)
525  c -- LOOP               END DO
526         yu10mx(j) = u10m(i,nsrf)            END DO
527         yu10my(j) = v10m(i,nsrf)  
        ywindsp(j) = SQRT(yu10mx(j)*yu10mx(j) + yu10my(j)*yu10my(j) )  
 c -- LOOP  
       END DO  
 C  
 C     IF bucket model for continent, copy soil water content  
       IF ( nsrf .eq. is_ter .and. .not. ok_veget ) THEN  
528            DO j = 1, knon            DO j = 1, knon
529              i = ni(j)               i = ni(j)
530              yqsol(j) = qsol(i)               DO k = 1, klev
531                    d_t(i, k) = d_t(i, k) + y_d_t(j, k)
532                    d_q(i, k) = d_q(i, k) + y_d_q(j, k)
533                    d_u(i, k) = d_u(i, k) + y_d_u(j, k)
534                    d_v(i, k) = d_v(i, k) + y_d_v(j, k)
535                    ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
536                 END DO
537            END DO            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  
538    
539        RETURN            ! diagnostic t, q a 2m et u, v a 10m
540        END  
541              DO j = 1, knon
542                 i = ni(j)
543                 uzon(j) = yu(j, 1) + y_d_u(j, 1)
544                 vmer(j) = yv(j, 1) + y_d_v(j, 1)
545                 tair1(j) = yt(j, 1) + y_d_t(j, 1)
546                 qair1(j) = yq(j, 1) + y_d_q(j, 1)
547                 zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
548                      1)))*(ypaprs(j, 1)-ypplay(j, 1))
549                 tairsol(j) = yts(j) + y_d_ts(j)
550                 rugo1(j) = yrugos(j)
551                 IF (nsrf == is_oce) THEN
552                    rugo1(j) = rugos(i, nsrf)
553                 END IF
554                 psfce(j) = ypaprs(j, 1)
555                 patm(j) = ypplay(j, 1)
556    
557                 qairsol(j) = yqsurf(j)
558              END DO
559    
560              CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &
561                   zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &
562                   yt10m, yq10m, yu10m, yustar)
563    
564              DO j = 1, knon
565                 i = ni(j)
566                 t2m(i, nsrf) = yt2m(j)
567                 q2m(i, nsrf) = yq2m(j)
568    
569                 ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
570                 u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
571                 v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)
572    
573              END DO
574    
575              CALL hbtm(knon, ypaprs, ypplay, yt2m, yq2m, yustar, y_flux_t, &
576                   y_flux_q, yu, yv, yt, yq, ypblh(:knon), ycapcl, yoliqcl, &
577                   ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
578    
579              DO j = 1, knon
580                 i = ni(j)
581                 pblh(i, nsrf) = ypblh(j)
582                 plcl(i, nsrf) = ylcl(j)
583                 capcl(i, nsrf) = ycapcl(j)
584                 oliqcl(i, nsrf) = yoliqcl(j)
585                 cteicl(i, nsrf) = ycteicl(j)
586                 pblt(i, nsrf) = ypblt(j)
587                 therm(i, nsrf) = ytherm(j)
588                 trmb1(i, nsrf) = ytrmb1(j)
589                 trmb2(i, nsrf) = ytrmb2(j)
590                 trmb3(i, nsrf) = ytrmb3(j)
591              END DO
592    
593              DO j = 1, knon
594                 DO k = 1, klev + 1
595                    i = ni(j)
596                    q2(i, k, nsrf) = yq2(j, k)
597                 END DO
598              END DO
599           end IF if_knon
600        END DO loop_surface
601    
602        ! On utilise les nouvelles surfaces
603    
604        rugos(:, is_oce) = rugmer
605        pctsrf = pctsrf_new
606    
607      END SUBROUTINE clmain
608    
609    end module clmain_m

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

  ViewVC Help
Powered by ViewVC 1.1.21