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

  ViewVC Help
Powered by ViewVC 1.1.21