/[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 305 by guez, Tue Sep 11 11:08:38 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_fall(klon) ! liquid water mass flux (kg / m2 / s), positive down
146        REAL cdragh(klon), cdragm(klon)      REAL ysnow_fall(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      yrugos = 0.
205        REAL y_d_ts(klon)      ypaprs = 0.
206        REAL y_d_t(klon, klev), y_d_q(klon, klev)      ypplay = 0.
207        REAL y_d_u(klon, klev), y_d_v(klon, klev)      ydelp = 0.
208        REAL y_flux_t(klon,klev), y_flux_q(klon,klev)      yrugoro = 0.
209        REAL y_flux_u(klon,klev), y_flux_v(klon,klev)      d_ts = 0.
210        REAL y_dflux_t(klon), y_dflux_q(klon)      flux_t = 0.
211        REAL ycoefh(klon,klev), ycoefm(klon,klev)      flux_q = 0.
212        REAL yu(klon,klev), yv(klon,klev)      flux_u = 0.
213        REAL yt(klon,klev), yq(klon,klev)      flux_v = 0.
214        REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)      fluxlat = 0.
215  c      d_t = 0.
216        LOGICAL ok_nonloc      d_q = 0.
217        PARAMETER (ok_nonloc=.FALSE.)      d_u = 0.
218        REAL ycoefm0(klon,klev), ycoefh0(klon,klev)      d_v = 0.
219        coefh = 0.
220  cIM 081204 hcl_Anne ? BEG      fqcalving = 0.
221        real yzlay(klon,klev),yzlev(klon,klev+1),yteta(klon,klev)      run_off_lic = 0.
222        real ykmm(klon,klev+1),ykmn(klon,klev+1)  
223        real ykmq(klon,klev+1)      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
224        real yq2(klon,klev+1),q2(klon,klev+1,nbsrf)      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
225        real q2diag(klon,klev+1)      ! (\`a affiner).
226  cIM 081204   real yustar(klon),y_cd_m(klon),y_cd_h(klon)  
227  cIM 081204 hcl_Anne ? END      pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
228  c      pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
229        REAL u1lay(klon), v1lay(klon)      pctsrf_pot(:, is_oce) = 1. - zmasq
230        REAL delp(klon,klev)      pctsrf_pot(:, is_sic) = 1. - zmasq
231        INTEGER i, k, nsrf  
232  cAA   INTEGER it      ! Tester si c'est le moment de lire le fichier:
233        INTEGER ni(klon), knon, j      if (mod(itap - 1, lmt_pas) == 0) then
234  c Introduction d'une variable "pourcentage potentiel" pour tenir compte         CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
235  c des eventuelles apparitions et/ou disparitions de la glace de mer      endif
236        REAL pctsrf_pot(klon,nbsrf)  
237              ! Boucler sur toutes les sous-fractions du sol:
238  c======================================================================  
239        REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.      loop_surface: DO nsrf = 1, nbsrf
240  c======================================================================         ! Define ni and knon:
241  c        
242  c maf pour sorties IOISPL en cas de debugagage         ni = 0
243  c         knon = 0
244        CHARACTER*80 cldebug  
245        SAVE cldebug         DO i = 1, klon
246        CHARACTER*8 cl_surf(nbsrf)            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
247        SAVE cl_surf            ! "potentielles"
248        INTEGER nhoridbg, nidbg            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
249        SAVE nhoridbg, nidbg               knon = knon + 1
250        INTEGER ndexbg(iim*(jjm+1))               ni(knon) = i
251        REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1), zjulian            END IF
252        REAL tabindx(klon)         END DO
253        REAL debugtab(iim,jjm+1)  
254        LOGICAL first_appel         if_knon: IF (knon /= 0) then
255        SAVE first_appel            DO j = 1, knon
256        DATA first_appel/.true./               i = ni(j)
257        LOGICAL:: debugindex = .false.               ypct(j) = pctsrf(i, nsrf)
258        integer idayref               yts(j) = ftsol(i, nsrf)
259        REAL t2m(klon,nbsrf), q2m(klon,nbsrf)               snow(j) = fsnow(i, nsrf)
260        REAL u10m(klon,nbsrf), v10m(klon,nbsrf)               yqsurf(j) = qsurf(i, nsrf)
261  c               yalb(j) = falbe(i, nsrf)
262        REAL yt2m(klon), yq2m(klon), yu10m(klon)               yrain_fall(j) = rain_fall(i)
263        REAL yustar(klon)               ysnow_fall(j) = snow_fall(i)
264  c -- LOOP               yagesno(j) = agesno(i, nsrf)
265         REAL yu10mx(klon)               yrugos(j) = frugs(i, nsrf)
266         REAL yu10my(klon)               yrugoro(j) = rugoro(i)
267         REAL ywindsp(klon)               yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)
268  c -- LOOP               ypaprs(j, klev + 1) = paprs(i, klev + 1)
269  c               y_run_off_lic_0(j) = run_off_lic_0(i)
270        REAL yt10m(klon), yq10m(klon)            END DO
271  cIM cf. AM : pbl, hbtm2 (Comme les autres diagnostics on cumule ds physic ce qui  
272  c   permet de sortir les grdeurs par sous surface)            ! For continent, copy soil water content
273        REAL pblh(klon,nbsrf)            IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
274        REAL plcl(klon,nbsrf)  
275        REAL capCL(klon,nbsrf)            ytsoil(:knon, :) = ftsoil(ni(:knon), :, nsrf)
276        REAL oliqCL(klon,nbsrf)  
277        REAL cteiCL(klon,nbsrf)            DO k = 1, klev
278        REAL pblT(klon,nbsrf)               DO j = 1, knon
279        REAL therm(klon,nbsrf)                  i = ni(j)
280        REAL trmb1(klon,nbsrf)                  ypaprs(j, k) = paprs(i, k)
281        REAL trmb2(klon,nbsrf)                  ypplay(j, k) = pplay(i, k)
282        REAL trmb3(klon,nbsrf)                  ydelp(j, k) = delp(i, k)
283        REAL ypblh(klon)                  yu(j, k) = u(i, k)
284        REAL ylcl(klon)                  yv(j, k) = v(i, k)
285        REAL ycapCL(klon)                  yt(j, k) = t(i, k)
286        REAL yoliqCL(klon)                  yq(j, k) = q(i, k)
287        REAL ycteiCL(klon)               END DO
       REAL ypblT(klon)  
       REAL ytherm(klon)  
       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)  
288            END DO            END DO
289            CALL histend(nidbg)  
290            CALL histsync(nidbg)            ! Calculer les géopotentiels de chaque couche:
291        ENDIF  
292              zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
293                   + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
294    
295              DO k = 2, klev
296                 zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
297                      * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
298                      * (ypplay(:knon, k - 1) - ypplay(:knon, k))
299              ENDDO
300    
301              CALL cdrag(nsrf, sqrt(yu(:knon, 1)**2 + yv(:knon, 1)**2), &
302                   yt(:knon, 1), yq(:knon, 1), zgeop(:knon, 1), ypaprs(:knon, 1), &
303                   yts(:knon), yqsurf(:knon), yrugos(:knon), ycdragm(:knon), &
304                   ycdragh(:knon))
305    
306              IF (iflag_pbl == 1) THEN
307                 ycdragm(:knon) = max(ycdragm(:knon), 0.)
308                 ycdragh(:knon) = max(ycdragh(:knon), 0.)
309              end IF
310    
311              ! on met un seuil pour ycdragm et ycdragh
312              IF (nsrf == is_oce) THEN
313                 ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
314                 ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
315              END IF
316    
317              IF (iflag_pbl >= 6) yq2(:knon, :) = q2(ni(:knon), :, nsrf)
318              call coef_diff_turb(nsrf, ni(:knon), ypaprs(:knon, :), &
319                   ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
320                   yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
321                   ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
322                        
323        DO k = 1, klev   ! epaisseur de couche            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
324        DO i = 1, klon                 ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
325           delp(i,k) = paprs(i,k)-paprs(i,k+1)                 ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
326        ENDDO                 y_flux_u(:knon))
327        ENDDO            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
328        DO i = 1, klon  ! vent de la premiere couche                 ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
329           zx_alf1 = 1.0                 ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
330           zx_alf2 = 1.0 - zx_alf1                 y_flux_v(:knon))
331           u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2  
332           v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2            CALL clqh(julien, nsrf, ni(:knon), ytsoil(:knon, :), yqsol(:knon), &
333        ENDDO                 mu0(ni(:knon)), yrugos(:knon), yrugoro(:knon), yu(:knon, 1), &
334  c                 yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), yt(:knon, :), &
335  c initialisation:                 yq(:knon, :), yts(:knon), ypaprs(:knon, :), ypplay(:knon, :), &
336  c                 ydelp(:knon, :), yrads(:knon), yalb(:knon), snow(:knon), &
337        DO i = 1, klon                 yqsurf(:knon), yrain_fall(:knon), ysnow_fall(:knon), &
338           rugmer(i) = 0.0                 yfluxlat(:knon), pctsrf_new_sic(ni(:knon)), yagesno(:knon), &
339           cdragh(i) = 0.0                 y_d_t(:knon, :), y_d_q(:knon, :), y_d_ts(:knon), &
340           cdragm(i) = 0.0                 yz0_new(:knon), y_flux_t(:knon), y_flux_q(:knon), &
341           dflux_t(i) = 0.0                 y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving(:knon), &
342           dflux_q(i) = 0.0                 y_ffonte(:knon), y_run_off_lic_0(:knon), y_run_off_lic(:knon))
343           zu1(i) = 0.0  
344           zv1(i) = 0.0            ! calculer la longueur de rugosite sur ocean
345        ENDDO  
346        ypct = 0.0            yrugm = 0.
347        yts = 0.0  
348        ysnow = 0.0            IF (nsrf == is_oce) THEN
349        yqsurf = 0.0               DO j = 1, knon
350        yalb = 0.0                  yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
351        yalblw = 0.0                       / rg + 0.11 * 14E-6 &
352        yrain_f = 0.0                       / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
353        ysnow_f = 0.0                  yrugm(j) = max(1.5E-05, yrugm(j))
354        yfder = 0.0               END DO
355        ytaux = 0.0            END IF
356        ytauy = 0.0  
357        ysolsw = 0.0            DO k = 1, klev
358        ysollw = 0.0               DO j = 1, knon
359        ysollwdown = 0.0                  i = ni(j)
360        yrugos = 0.0                  y_d_t(j, k) = y_d_t(j, k) * ypct(j)
361        yu1 = 0.0                  y_d_q(j, k) = y_d_q(j, k) * ypct(j)
362        yv1 = 0.0                  y_d_u(j, k) = y_d_u(j, k) * ypct(j)
363        yrads = 0.0                  y_d_v(j, k) = y_d_v(j, k) * ypct(j)
364        ypaprs = 0.0               END DO
365        ypplay = 0.0            END DO
366        ydelp = 0.0  
367        yu = 0.0            flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
368        yv = 0.0            flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
369        yt = 0.0            flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
370        yq = 0.0            flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
371        pctsrf_new = 0.0  
372        y_flux_u = 0.0            falbe(:, nsrf) = 0.
373        y_flux_v = 0.0            fsnow(:, nsrf) = 0.
374  C$$ PB            qsurf(:, nsrf) = 0.
375        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  
376            DO j = 1, knon            DO j = 1, knon
377              i = ni(j)               i = ni(j)
378              yqsol(j) = qsol(i)               d_ts(i, nsrf) = y_d_ts(j)
379                 falbe(i, nsrf) = yalb(j)
380                 fsnow(i, nsrf) = snow(j)
381                 qsurf(i, nsrf) = yqsurf(j)
382                 frugs(i, nsrf) = yz0_new(j)
383                 fluxlat(i, nsrf) = yfluxlat(j)
384                 IF (nsrf == is_oce) THEN
385                    rugmer(i) = yrugm(j)
386                    frugs(i, nsrf) = yrugm(j)
387                 END IF
388                 agesno(i, nsrf) = yagesno(j)
389                 fqcalving(i, nsrf) = y_fqcalving(j)
390                 ffonte(i, nsrf) = y_ffonte(j)
391                 cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
392                 cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
393                 dflux_t(i) = dflux_t(i) + y_dflux_t(j) * ypct(j)
394                 dflux_q(i) = dflux_q(i) + y_dflux_q(j) * ypct(j)
395            END DO            END DO
396        ELSE            IF (nsrf == is_ter) THEN
397            yqsol(:)=0.               qsol(ni(:knon)) = yqsol(:knon)
398        ENDIF            else IF (nsrf == is_lic) THEN
399  c$$$ PB ajour pour soil               DO j = 1, knon
400        DO k = 1, nsoilmx                  i = ni(j)
401          DO j = 1, knon                  run_off_lic_0(i) = y_run_off_lic_0(j)
402            i = ni(j)                  run_off_lic(i) = y_run_off_lic(j)
403            ytsoil(j,k) = ftsoil(i,k,nsrf)               END DO
404          END DO              END IF
405        END DO  
406        DO k = 1, klev            ftsoil(:, :, nsrf) = 0.
407        DO j = 1, knon            ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)
408        i = ni(j)  
409          ypaprs(j,k) = paprs(i,k)            DO j = 1, knon
410          ypplay(j,k) = pplay(i,k)               i = ni(j)
411          ydelp(j,k) = delp(i,k)               DO k = 1, klev
412          yu(j,k) = u(i,k)                  d_t(i, k) = d_t(i, k) + y_d_t(j, k)
413          yv(j,k) = v(i,k)                  d_q(i, k) = d_q(i, k) + y_d_q(j, k)
414          yt(j,k) = t(i,k)                  d_u(i, k) = d_u(i, k) + y_d_u(j, k)
415          yq(j,k) = q(i,k)                  d_v(i, k) = d_v(i, k) + y_d_v(j, k)
416        ENDDO               END DO
417        ENDDO            END DO
418  c  
419  c            forall (k = 2:klev) coefh(ni(:knon), k) &
420  c calculer Cdrag et les coefficients d'echange                 = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypct(:knon)
421        CALL coefkz(nsrf, knon, ypaprs, ypplay,  
422  cIM 261103            ! diagnostic t, q a 2m et u, v a 10m
423       .     ksta, ksta_ter,  
424  cIM 261103            DO j = 1, knon
425       .            yts, yrugos, yu, yv, yt, yq,               i = ni(j)
426       .            yqsurf,               u1(j) = yu(j, 1) + y_d_u(j, 1)
427       .            ycoefm, ycoefh)               v1(j) = yv(j, 1) + y_d_v(j, 1)
428  cIM 081204 BEG               tair1(j) = yt(j, 1) + y_d_t(j, 1)
429  cCR test               qair1(j) = yq(j, 1) + y_d_q(j, 1)
430        if (iflag_pbl.eq.1) then               zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
431  cIM 081204 END                    1))) * (ypaprs(j, 1)-ypplay(j, 1))
432          CALL coefkz2(nsrf, knon, ypaprs, ypplay,yt,               tairsol(j) = yts(j) + y_d_ts(j)
433       .                  ycoefm0, ycoefh0)               rugo1(j) = yrugos(j)
434          DO k = 1, klev               IF (nsrf == is_oce) THEN
435          DO i = 1, knon                  rugo1(j) = frugs(i, nsrf)
436             ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))               END IF
437             ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))               psfce(j) = ypaprs(j, 1)
438          ENDDO               patm(j) = ypplay(j, 1)
439          ENDDO            END DO
440        endif  
441  c            CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), qair1, &
442  cIM cf JLD : on seuille ycoefm et ycoefh                 zgeo1, tairsol, yqsurf(:knon), rugo1, psfce, patm, yt2m, yq2m, &
443        if (nsrf.eq.is_oce) then                 yt10m, yq10m, wind10m(:knon), ustar(:knon))
444           do j=1,knon  
445  c           ycoefm(j,1)=min(ycoefm(j,1),1.1E-3)            DO j = 1, knon
446              ycoefm(j,1)=min(ycoefm(j,1),cdmmax)               i = ni(j)
447  c           ycoefh(j,1)=min(ycoefh(j,1),1.1E-3)               t2m(i, nsrf) = yt2m(j)
448              ycoefh(j,1)=min(ycoefh(j,1),cdhmax)               q2m(i, nsrf) = yq2m(j)
449           enddo  
450        endif               u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
451                      / sqrt(u1(j)**2 + v1(j)**2)
452  c               v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
453  cIM: 261103                    / sqrt(u1(j)**2 + v1(j)**2)
454        if (ok_kzmin) THEN            END DO
455  cIM cf FH: 201103 BEG  
456  c   Calcul d'une diffusion minimale pour les conditions tres stables.            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
457        call coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycoefm                 y_flux_q(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
458       .   ,ycoefm0,ycoefh0)                 yq(:knon, :), ypblh(:knon), ycapcl, yoliqcl, ycteicl, ypblt, &
459  c      call dump2d(iim,jjm-1,ycoefm(2:klon-1,2), 'KZ         ')                 ytherm, ylcl)
460  c      call dump2d(iim,jjm-1,ycoefm0(2:klon-1,2),'KZMIN      ')  
461              DO j = 1, knon
462         if ( 1.eq.1 ) then               i = ni(j)
463         DO k = 1, klev               pblh(i, nsrf) = ypblh(j)
464         DO i = 1, knon               plcl(i, nsrf) = ylcl(j)
465            ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))               capcl(i, nsrf) = ycapcl(j)
466            ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))               oliqcl(i, nsrf) = yoliqcl(j)
467         ENDDO               cteicl(i, nsrf) = ycteicl(j)
468         ENDDO               pblt(i, nsrf) = ypblt(j)
469         endif               therm(i, nsrf) = ytherm(j)
470  cIM cf FH: 201103 END            END DO
471        endif !ok_kzmin  
472  cIM: 261103            IF (iflag_pbl >= 6) q2(ni(:knon), :, nsrf) = yq2(:knon, :)
473           else
474              fsnow(:, nsrf) = 0.
475        IF (iflag_pbl.ge.3) then         end IF if_knon
476        END DO loop_surface
477  cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
478  c MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin      ! On utilise les nouvelles surfaces
479  cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      frugs(:, is_oce) = rugmer
480        pctsrf(:, is_oce) = pctsrf_new_oce
481           yzlay(1:knon,1)=      pctsrf(:, is_sic) = pctsrf_new_sic
482       .   RD*yt(1:knon,1)/(0.5*(ypaprs(1:knon,1)+ypplay(1:knon,1)))  
483       .   *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG      CALL histwrite_phy("run_off_lic", run_off_lic)
484           do k=2,klev  
485              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  
486    
487        RETURN  end module pbl_surface_m
       END  

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

  ViewVC Help
Powered by ViewVC 1.1.21