/[lmdze]/trunk/phylmd/Interface_surf/pbl_surface.f
ViewVC logotype

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

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

  ViewVC Help
Powered by ViewVC 1.1.21