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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21