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

Legend:
Removed from v.13  
changed lines
  Added in v.175

  ViewVC Help
Powered by ViewVC 1.1.21