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

Legend:
Removed from v.14  
changed lines
  Added in v.237

  ViewVC Help
Powered by ViewVC 1.1.21