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

Legend:
Removed from v.7  
changed lines
  Added in v.227

  ViewVC Help
Powered by ViewVC 1.1.21