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

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

  ViewVC Help
Powered by ViewVC 1.1.21