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

Diff of /trunk/phylmd/Interface_surf/pbl_surface.f90

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/Sources/phylmd/clmain.f revision 225 by guez, Mon Oct 16 12:35:41 2017 UTC trunk/phylmd/Interface_surf/pbl_surface.f revision 324 by guez, Wed Feb 6 15:58:03 2019 UTC
# Line 1  Line 1 
1  module clmain_m  module pbl_surface_m
2    
3    IMPLICIT NONE    IMPLICIT NONE
4    
5  contains  contains
6    
7    SUBROUTINE clmain(dtime, pctsrf, t, q, u, v, julien, mu0, ftsol, cdmmax, &    SUBROUTINE pbl_surface(pctsrf, t, q, u, v, julien, mu0, ftsol, cdmmax, &
8         cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, pplay, fsnow, &         cdhmax, ftsoil, qsol, paprs, play, fsnow, fqsurf, falbe, fluxlat, &
9         qsurf, evap, falbe, fluxlat, rain_fall, snow_f, fsolsw, fsollw, frugs, &         rain_fall, snow_fall, frugs, agesno, rugoro, d_t, d_q, d_u, d_v, d_ts, &
10         agesno, rugoro, d_t, d_q, d_u, d_v, d_ts, flux_t, flux_q, flux_u, &         flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, dflux_t, dflux_q, &
11         flux_v, cdragh, cdragm, q2, dflux_t, dflux_q, ycoefh, zu1, zv1, t2m, &         coefh, t2m, q2m, u10m_srf, v10m_srf, pblh, capcl, oliqcl, cteicl, pblt, &
12         q2m, u10m_srf, v10m_srf, pblh, capcl, oliqcl, cteicl, pblt, therm, &         therm, plcl, fqcalving, ffonte, run_off_lic_0, albsol, sollw, solsw, &
13         trmb1, trmb2, trmb3, plcl, fqcalving, ffonte, run_off_lic_0)         tsol)
14    
15      ! From phylmd/clmain.F, version 1.6, 2005/11/16 14:47:19      ! From phylmd/clmain.F, version 1.6, 2005/11/16 14:47:19
16      ! Author: Z. X. Li (LMD/CNRS), date: 1993/08/18      ! Author: Z. X. Li (LMD/CNRS)
17        ! Date: Aug. 18th, 1993
18      ! Objet : interface de couche limite (diffusion verticale)      ! Objet : interface de couche limite (diffusion verticale)
19    
20      ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul      ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
# Line 21  contains Line 22  contains
22      ! ne tient pas compte de la diff\'erentiation des sous-fractions      ! ne tient pas compte de la diff\'erentiation des sous-fractions
23      ! de sol.      ! de sol.
24    
25      ! Pour pouvoir extraire les coefficients d'\'echanges et le vent      use cdrag_m, only: cdrag
     ! dans la premi\`ere couche, trois champs ont \'et\'e cr\'e\'es : "ycoefh",  
     ! "zu1" et "zv1". Nous avons moyenn\'e les valeurs de ces trois  
     ! champs sur les quatre sous-surfaces du mod\`ele.  
   
26      use clqh_m, only: clqh      use clqh_m, only: clqh
27      use clvent_m, only: clvent      use clvent_m, only: clvent
28      use coefkz_m, only: coefkz      use coef_diff_turb_m, only: coef_diff_turb
29      use coefkzmin_m, only: coefkzmin      USE conf_gcm_m, ONLY: lmt_pas
     USE conf_gcm_m, ONLY: prt_level, lmt_pas  
30      USE conf_phys_m, ONLY: iflag_pbl      USE conf_phys_m, ONLY: iflag_pbl
31      USE dimphy, ONLY: klev, klon, zmasq      USE dimphy, ONLY: klev, klon
32      USE dimsoil, ONLY: nsoilmx      USE dimsoil, ONLY: nsoilmx
33      use hbtm_m, only: hbtm      use hbtm_m, only: hbtm
34        USE histwrite_phy_m, ONLY: histwrite_phy
35      USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf      USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
36      USE interfoce_lim_m, ONLY: interfoce_lim      USE interfoce_lim_m, ONLY: interfoce_lim
37        use phyetat0_m, only: masque
38      use stdlevvar_m, only: stdlevvar      use stdlevvar_m, only: stdlevvar
39      USE suphec_m, ONLY: rd, rg, rkappa      USE suphec_m, ONLY: rd, rg, rsigma
40      use time_phylmdz, only: itap      use time_phylmdz, only: itap
     use ustarhb_m, only: ustarhb  
     use vdif_kcay_m, only: vdif_kcay  
     use yamada4_m, only: yamada4  
   
     REAL, INTENT(IN):: dtime ! interval du temps (secondes)  
41    
42      REAL, INTENT(inout):: pctsrf(klon, nbsrf)      REAL, INTENT(inout):: pctsrf(klon, nbsrf)
43      ! tableau des pourcentages de surface de chaque maille      ! pourcentages de surface de chaque maille
44    
45      REAL, INTENT(IN):: t(klon, klev) ! temperature (K)      REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
46      REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg / kg)      REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg / kg)
47      REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse      REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
48      INTEGER, INTENT(IN):: julien ! jour de l'annee en cours      INTEGER, INTENT(IN):: julien ! jour de l'annee en cours
49      REAL, intent(in):: mu0(klon) ! cosinus de l'angle solaire zenithal          REAL, intent(in):: mu0(klon) ! cosinus de l'angle solaire zenithal    
50      REAL, INTENT(IN):: ftsol(:, :) ! (klon, nbsrf) temp\'erature du sol (en K)  
51        REAL, INTENT(IN):: ftsol(:, :) ! (klon, nbsrf)
52        ! skin temperature of surface fraction, in K
53    
54      REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh      REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
     REAL, INTENT(IN):: ksta, ksta_ter  
     LOGICAL, INTENT(IN):: ok_kzmin  
55    
56      REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)      REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
57      ! soil temperature of surface fraction      ! soil temperature of surface fraction
# Line 66  contains Line 60  contains
60      ! column-density of water in soil, in kg m-2      ! column-density of water in soil, in kg m-2
61    
62      REAL, INTENT(IN):: paprs(klon, klev + 1) ! pression a intercouche (Pa)      REAL, INTENT(IN):: paprs(klon, klev + 1) ! pression a intercouche (Pa)
63      REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)      REAL, INTENT(IN):: play(klon, klev) ! pression au milieu de couche (Pa)
64      REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse      REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse
65      REAL qsurf(klon, nbsrf)      REAL, INTENT(inout):: fqsurf(klon, nbsrf)
     REAL evap(klon, nbsrf)  
66      REAL, intent(inout):: falbe(klon, nbsrf)      REAL, intent(inout):: falbe(klon, nbsrf)
67    
68      REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)      REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)
69        ! flux de chaleur latente, en W m-2
70    
71      REAL, intent(in):: rain_fall(klon)      REAL, intent(in):: rain_fall(klon)
72      ! liquid water mass flux (kg / m2 / s), positive down      ! liquid water mass flux (kg / m2 / s), positive down
73    
74      REAL, intent(in):: snow_f(klon)      REAL, intent(in):: snow_fall(klon)
75      ! solid water mass flux (kg / m2 / s), positive down      ! solid water mass flux (kg / m2 / s), positive down
76    
     REAL, INTENT(IN):: fsolsw(klon, nbsrf), fsollw(klon, nbsrf)  
77      REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)      REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)
78      real agesno(klon, nbsrf)      real agesno(klon, nbsrf)
79      REAL, INTENT(IN):: rugoro(klon)      REAL, INTENT(IN):: rugoro(klon)
80    
81      REAL d_t(klon, klev), d_q(klon, klev)      REAL, intent(out):: d_t(:, :), d_q(:, :) ! (klon, klev)
82      ! d_t------output-R- le changement pour "t"      ! changement pour t et q
     ! d_q------output-R- le changement pour "q"  
83    
84      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
85      ! changement pour "u" et "v"      ! changement pour "u" et "v"
# Line 94  contains Line 87  contains
87      REAL, intent(out):: d_ts(:, :) ! (klon, nbsrf) variation of ftsol      REAL, intent(out):: d_ts(:, :) ! (klon, nbsrf) variation of ftsol
88    
89      REAL, intent(out):: flux_t(klon, nbsrf)      REAL, intent(out):: flux_t(klon, nbsrf)
90      ! flux de chaleur sensible (Cp T) (W / m2) (orientation positive vers      ! flux de chaleur sensible (c_p T) (W / m2) (orientation positive
91      ! le bas) à la surface      ! vers le bas) à la surface
92    
93      REAL, intent(out):: flux_q(klon, nbsrf)      REAL, intent(out):: flux_q(klon, nbsrf)
94      ! flux de vapeur d'eau (kg / m2 / s) à la surface      ! flux de vapeur d'eau (kg / m2 / s) à la surface
95    
96      REAL, intent(out):: flux_u(klon, nbsrf), flux_v(klon, nbsrf)      REAL, intent(out):: flux_u(:, :), flux_v(:, :) ! (klon, nbsrf)
97      ! tension du vent à la surface, en Pa      ! tension du vent (flux turbulent de vent) à la surface, en Pa
98    
99      REAL, INTENT(out):: cdragh(klon), cdragm(klon)      REAL, INTENT(out):: cdragh(klon), cdragm(klon)
100      real q2(klon, klev + 1, nbsrf)      real q2(klon, klev + 1, nbsrf)
101    
102      REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)      ! Ocean slab:
103      ! dflux_t derive du flux sensible      REAL, INTENT(out):: dflux_t(klon) ! derive du flux sensible
104      ! dflux_q derive du flux latent      REAL, INTENT(out):: dflux_q(klon) ! derive du flux latent
105      ! IM "slab" ocean  
106        REAL, intent(out):: coefh(:, 2:) ! (klon, 2:klev)
107        ! Pour pouvoir extraire les coefficients d'\'echange, le champ
108        ! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
109        ! ce champ sur les quatre sous-surfaces du mod\`ele.
110    
     REAL, intent(out):: ycoefh(klon, klev)  
     REAL, intent(out):: zu1(klon), zv1(klon)  
111      REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)      REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)
112    
113      REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)      REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)
# Line 127  contains Line 122  contains
122      REAL cteicl(klon, nbsrf)      REAL cteicl(klon, nbsrf)
123      REAL, INTENT(inout):: pblt(klon, nbsrf) ! T au nveau HCL      REAL, INTENT(inout):: pblt(klon, nbsrf) ! T au nveau HCL
124      REAL therm(klon, nbsrf)      REAL therm(klon, nbsrf)
     REAL trmb1(klon, nbsrf)  
     ! trmb1-------deep_cape  
     REAL trmb2(klon, nbsrf)  
     ! trmb2--------inhibition  
     REAL trmb3(klon, nbsrf)  
     ! trmb3-------Point Omega  
125      REAL plcl(klon, nbsrf)      REAL plcl(klon, nbsrf)
126      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)  
127      ! ffonte----Flux thermique utilise pour fondre la neige      REAL, intent(out):: fqcalving(klon, nbsrf)
128      ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la      ! flux d'eau "perdue" par la surface et necessaire pour limiter la
129      !           hauteur de neige, en kg / m2 / s      ! hauteur de neige, en kg / m2 / s
130      REAL run_off_lic_0(klon)  
131        real ffonte(klon, nbsrf) ! flux thermique utilise pour fondre la neige
132        REAL, intent(inout):: run_off_lic_0(:) ! (klon)
133    
134        REAL, intent(out):: albsol(:) ! (klon)
135        ! albedo du sol total, visible, moyen par maille
136    
137        REAL, intent(in):: sollw(:) ! (klon)
138        ! surface net downward longwave flux, in W m-2
139    
140        REAL, intent(in):: solsw(:) ! (klon)
141        ! surface net downward shortwave flux, in W m-2
142    
143        REAL, intent(in):: tsol(:) ! (klon)
144    
145      ! Local:      ! Local:
146    
147      LOGICAL:: firstcal = .true.      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
148        REAL fsolsw(klon, nbsrf) ! flux solaire absorb\'e pour chaque sous-surface
149    
150      ! la nouvelle repartition des surfaces sortie de l'interface      ! la nouvelle repartition des surfaces sortie de l'interface
151      REAL, save:: pctsrf_new_oce(klon)      REAL, save:: pctsrf_new_oce(klon)
152      REAL, save:: pctsrf_new_sic(klon)      REAL, save:: pctsrf_new_sic(klon)
153    
154      REAL y_fqcalving(klon), y_ffonte(klon)      REAL y_fqcalving(klon), y_ffonte(klon)
155      real y_run_off_lic_0(klon)      real y_run_off_lic_0(klon), y_run_off_lic(klon)
156        REAL run_off_lic(klon) ! ruissellement total
157      REAL rugmer(klon)      REAL rugmer(klon)
158      REAL ytsoil(klon, nsoilmx)      REAL ytsoil(klon, nsoilmx)
159      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      REAL yts(klon), ypctsrf(klon), yz0_new(klon)
160        real yrugos(klon) ! longueur de rugosite (en m)
161      REAL yalb(klon)      REAL yalb(klon)
   
     REAL u1lay(klon), v1lay(klon) ! vent dans la premi\`ere couche, pour  
                               ! une sous-surface donnée  
       
162      REAL snow(klon), yqsurf(klon), yagesno(klon)      REAL snow(klon), yqsurf(klon), yagesno(klon)
163      real yqsol(klon) ! column-density of water in soil, in kg m-2      real yqsol(klon) ! column-density of water in soil, in kg m-2
164      REAL yrain_f(klon) ! liquid water mass flux (kg / m2 / s), positive down      REAL yrain_fall(klon) ! liquid water mass flux (kg / m2 / s), positive down
165      REAL ysnow_f(klon) ! solid water mass flux (kg / m2 / s), positive down      REAL ysnow_fall(klon) ! solid water mass flux (kg / m2 / s), positive down
166      REAL yrugm(klon), yrads(klon), yrugoro(klon)      REAL yrugm(klon), radsol(klon), yrugoro(klon)
167      REAL yfluxlat(klon)      REAL yfluxlat(klon)
168      REAL y_d_ts(klon)      REAL y_d_ts(klon)
169      REAL y_d_t(klon, klev), y_d_q(klon, klev)      REAL y_d_t(klon, klev), y_d_q(klon, klev)
# Line 170  contains Line 171  contains
171      REAL y_flux_t(klon), y_flux_q(klon)      REAL y_flux_t(klon), y_flux_q(klon)
172      REAL y_flux_u(klon), y_flux_v(klon)      REAL y_flux_u(klon), y_flux_v(klon)
173      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
174      REAL coefh(klon, klev), coefm(klon, klev)      REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
175        real ycdragh(klon), ycdragm(klon)
176      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
177      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
178      REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)      REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
   
     REAL ycoefm0(klon, klev), ycoefh0(klon, klev)  
   
     REAL yzlay(klon, klev), yzlev(klon, klev + 1), yteta(klon, klev)  
     REAL ykmm(klon, klev + 1), ykmn(klon, klev + 1)  
     REAL ykmq(klon, klev + 1)  
179      REAL yq2(klon, klev + 1)      REAL yq2(klon, klev + 1)
     REAL q2diag(klon, klev + 1)  
   
180      REAL delp(klon, klev)      REAL delp(klon, klev)
181      INTEGER i, k, nsrf      INTEGER i, k, nsrf
   
182      INTEGER ni(klon), knon, j      INTEGER ni(klon), knon, j
183    
184      REAL pctsrf_pot(klon, nbsrf)      REAL pctsrf_pot(klon, nbsrf)
185      ! "pourcentage potentiel" pour tenir compte des \'eventuelles      ! "pourcentage potentiel" pour tenir compte des \'eventuelles
186      ! apparitions ou disparitions de la glace de mer      ! apparitions ou disparitions de la glace de mer
187    
188      REAL yt2m(klon), yq2m(klon), yu10m(klon)      REAL yt2m(klon), yq2m(klon), wind10m(klon)
189      REAL yustar(klon)      REAL ustar(klon)
190    
191      REAL yt10m(klon), yq10m(klon)      REAL yt10m(klon), yq10m(klon)
192      REAL ypblh(klon)      REAL ypblh(klon)
# Line 203  contains Line 196  contains
196      REAL ycteicl(klon)      REAL ycteicl(klon)
197      REAL ypblt(klon)      REAL ypblt(klon)
198      REAL ytherm(klon)      REAL ytherm(klon)
199      REAL ytrmb1(klon)      REAL u1(klon), v1(klon)
     REAL ytrmb2(klon)  
     REAL ytrmb3(klon)  
     REAL uzon(klon), vmer(klon)  
200      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
201      REAL psfce(klon), patm(klon)      REAL psfce(klon), patm(klon)
202        REAL zgeo1(klon)
     REAL qairsol(klon), zgeo1(klon)  
203      REAL rugo1(klon)      REAL rugo1(klon)
204        REAL zgeop(klon, klev)
     ! utiliser un jeu de fonctions simples                
     LOGICAL zxli  
     PARAMETER (zxli=.FALSE.)  
205    
206      !------------------------------------------------------------      !------------------------------------------------------------
207    
208        albsol = sum(falbe * pctsrf, dim = 2)
209    
210        ! R\'epartition sous maille des flux longwave et shortwave
211        ! R\'epartition du longwave par sous-surface lin\'earis\'ee
212    
213        forall (nsrf = 1:nbsrf)
214           fsollw(:, nsrf) = sollw + 4. * RSIGMA * tsol**3 &
215                * (tsol - ftsol(:, nsrf))
216           fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
217        END forall
218    
219      ytherm = 0.      ytherm = 0.
220    
221      DO k = 1, klev ! epaisseur de couche      DO k = 1, klev ! epaisseur de couche
# Line 233  contains Line 230  contains
230      cdragm = 0.      cdragm = 0.
231      dflux_t = 0.      dflux_t = 0.
232      dflux_q = 0.      dflux_q = 0.
     zu1 = 0.  
     zv1 = 0.  
     ypct = 0.  
     yqsurf = 0.  
     yrain_f = 0.  
     ysnow_f = 0.  
233      yrugos = 0.      yrugos = 0.
234      ypaprs = 0.      ypaprs = 0.
235      ypplay = 0.      ypplay = 0.
236      ydelp = 0.      ydelp = 0.
     yu = 0.  
     yv = 0.  
     yt = 0.  
     yq = 0.  
     y_dflux_t = 0.  
     y_dflux_q = 0.  
237      yrugoro = 0.      yrugoro = 0.
238      d_ts = 0.      d_ts = 0.
239      flux_t = 0.      flux_t = 0.
# Line 260  contains Line 245  contains
245      d_q = 0.      d_q = 0.
246      d_u = 0.      d_u = 0.
247      d_v = 0.      d_v = 0.
248      ycoefh = 0.      coefh = 0.
249        fqcalving = 0.
250        run_off_lic = 0.
251    
252      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
253      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
254      ! (\`a affiner)      ! (\`a affiner).
255    
256      pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)      pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
257      pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)      pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
258      pctsrf_pot(:, is_oce) = 1. - zmasq      pctsrf_pot(:, is_oce) = 1. - masque
259      pctsrf_pot(:, is_sic) = 1. - zmasq      pctsrf_pot(:, is_sic) = 1. - masque
260    
261      ! Tester si c'est le moment de lire le fichier:      ! Tester si c'est le moment de lire le fichier:
262      if (mod(itap - 1, lmt_pas) == 0) then      if (mod(itap - 1, lmt_pas) == 0) then
# Line 279  contains Line 266  contains
266      ! Boucler sur toutes les sous-fractions du sol:      ! Boucler sur toutes les sous-fractions du sol:
267    
268      loop_surface: DO nsrf = 1, nbsrf      loop_surface: DO nsrf = 1, nbsrf
269         ! Chercher les indices :         ! Define ni and knon:
270    
271         ni = 0         ni = 0
272         knon = 0         knon = 0
273    
274         DO i = 1, klon         DO i = 1, klon
275            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
276            ! "potentielles"            ! "potentielles"
# Line 292  contains Line 281  contains
281         END DO         END DO
282    
283         if_knon: IF (knon /= 0) then         if_knon: IF (knon /= 0) then
284            DO j = 1, knon            ypctsrf(:knon) = pctsrf(ni(:knon), nsrf)
285               i = ni(j)            yts(:knon) = ftsol(ni(:knon), nsrf)
286               ypct(j) = pctsrf(i, nsrf)            snow(:knon) = fsnow(ni(:knon), nsrf)
287               yts(j) = ftsol(i, nsrf)            yqsurf(:knon) = fqsurf(ni(:knon), nsrf)
288               snow(j) = fsnow(i, nsrf)            yalb(:knon) = falbe(ni(:knon), nsrf)
289               yqsurf(j) = qsurf(i, nsrf)            yrain_fall(:knon) = rain_fall(ni(:knon))
290               yalb(j) = falbe(i, nsrf)            ysnow_fall(:knon) = snow_fall(ni(:knon))
291               yrain_f(j) = rain_fall(i)            yagesno(:knon) = agesno(ni(:knon), nsrf)
292               ysnow_f(j) = snow_f(i)            yrugos(:knon) = frugs(ni(:knon), nsrf)
293               yagesno(j) = agesno(i, nsrf)            yrugoro(:knon) = rugoro(ni(:knon))
294               yrugos(j) = frugs(i, nsrf)            radsol(:knon) = fsolsw(ni(:knon), nsrf) + fsollw(ni(:knon), nsrf)
295               yrugoro(j) = rugoro(i)            ypaprs(:knon, klev + 1) = paprs(ni(:knon), klev + 1)
296               u1lay(j) = u(i, 1)            y_run_off_lic_0(:knon) = run_off_lic_0(ni(:knon))
              v1lay(j) = v(i, 1)  
              yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)  
              ypaprs(j, klev + 1) = paprs(i, klev + 1)  
              y_run_off_lic_0(j) = run_off_lic_0(i)  
           END DO  
297    
298            ! For continent, copy soil water content            ! For continent, copy soil water content
299            IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))            IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
# Line 320  contains Line 304  contains
304               DO j = 1, knon               DO j = 1, knon
305                  i = ni(j)                  i = ni(j)
306                  ypaprs(j, k) = paprs(i, k)                  ypaprs(j, k) = paprs(i, k)
307                  ypplay(j, k) = pplay(i, k)                  ypplay(j, k) = play(i, k)
308                  ydelp(j, k) = delp(i, k)                  ydelp(j, k) = delp(i, k)
309                  yu(j, k) = u(i, k)                  yu(j, k) = u(i, k)
310                  yv(j, k) = v(i, k)                  yv(j, k) = v(i, k)
# Line 329  contains Line 313  contains
313               END DO               END DO
314            END DO            END DO
315    
316            ! calculer Cdrag et les coefficients d'echange            ! Calculer les géopotentiels de chaque couche:
           CALL coefkz(nsrf, ypaprs, ypplay, ksta, ksta_ter, yts(:knon), &  
                yrugos, yu, yv, yt, yq, yqsurf(:knon), coefm(:knon, :), &  
                coefh(:knon, :))  
           IF (iflag_pbl == 1) THEN  
              CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)  
              coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))  
              coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))  
           END IF  
317    
318            ! on met un seuil pour coefm et coefh            zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
319            IF (nsrf == is_oce) THEN                 + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
              coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)  
              coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)  
           END IF  
320    
321            IF (ok_kzmin) THEN            DO k = 2, klev
322               ! Calcul d'une diffusion minimale pour les conditions tres stables               zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
323               CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &                    * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
324                    coefm(:knon, 1), ycoefm0, ycoefh0)                    * (ypplay(:knon, k - 1) - ypplay(:knon, k))
325               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))            ENDDO
326               coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))  
327            END IF            CALL cdrag(nsrf, sqrt(yu(:knon, 1)**2 + yv(:knon, 1)**2), &
328                   yt(:knon, 1), yq(:knon, 1), zgeop(:knon, 1), ypaprs(:knon, 1), &
329                   yts(:knon), yqsurf(:knon), yrugos(:knon), ycdragm(:knon), &
330                   ycdragh(:knon))
331    
332            IF (iflag_pbl >= 3) THEN            IF (iflag_pbl == 1) THEN
333               ! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et               ycdragm(:knon) = max(ycdragm(:knon), 0.)
334               ! Fr\'ed\'eric Hourdin               ycdragh(:knon) = max(ycdragh(:knon), 0.)
335               yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &            end IF
                   + ypplay(:knon, 1))) &  
                   * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg  
              DO k = 2, klev  
                 yzlay(1:knon, k) = yzlay(1:knon, k-1) &  
                      + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &  
                      / ypaprs(1:knon, k) &  
                      * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg  
              END DO  
              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))  
              END DO  
              yzlev(1:knon, 1) = 0.  
              yzlev(:knon, klev + 1) = 2. * yzlay(:knon, klev) &  
                   - yzlay(:knon, klev - 1)  
              DO k = 2, klev  
                 yzlev(1:knon, k) = 0.5 * (yzlay(1:knon, k) + yzlay(1:knon, k-1))  
              END DO  
              DO k = 1, klev + 1  
                 DO j = 1, knon  
                    i = ni(j)  
                    yq2(j, k) = q2(i, k, nsrf)  
                 END DO  
              END DO  
   
              CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)  
              IF (prt_level > 9) PRINT *, 'USTAR = ', yustar  
   
              ! iflag_pbl peut \^etre utilis\'e comme longueur de m\'elange  
   
              IF (iflag_pbl >= 11) THEN  
                 CALL vdif_kcay(knon, dtime, rg, ypaprs, yzlev, yzlay, yu, yv, &  
                      yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, yustar, &  
                      iflag_pbl)  
              ELSE  
                 CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &  
                      coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)  
              END IF  
336    
337               coefm(:knon, 2:) = ykmm(:knon, 2:klev)            ! on met un seuil pour ycdragm et ycdragh
338               coefh(:knon, 2:) = ykmn(:knon, 2:klev)            IF (nsrf == is_oce) THEN
339                 ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
340                 ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
341            END IF            END IF
342    
343            ! calculer la diffusion des vitesses "u" et "v"            IF (iflag_pbl >= 6) yq2(:knon, :) = q2(ni(:knon), :, nsrf)
344            CALL clvent(knon, dtime, u1lay(:knon), v1lay(:knon), &            call coef_diff_turb(nsrf, ni(:knon), ypaprs(:knon, :), &
345                 coefm(:knon, :), yt, yu, ypaprs, ypplay, ydelp, y_d_u, &                 ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
346                   yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
347                   ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
348    
349              CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
350                   ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
351                   ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
352                 y_flux_u(:knon))                 y_flux_u(:knon))
353            CALL clvent(knon, dtime, u1lay(:knon), v1lay(:knon), &            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
354                 coefm(:knon, :), yt, yv, ypaprs, ypplay, ydelp, y_d_v, &                 ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
355                   ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
356                 y_flux_v(:knon))                 y_flux_v(:knon))
357    
358            ! calculer la diffusion de "q" et de "h"            CALL clqh(julien, nsrf, ni(:knon), ytsoil(:knon, :), yqsol(:knon), &
359            CALL clqh(dtime, julien, firstcal, nsrf, ni(:knon), &                 mu0(ni(:knon)), yrugos(:knon), yrugoro(:knon), yu(:knon, 1), &
360                 ytsoil(:knon, :), yqsol(:knon), mu0, yrugos, yrugoro, &                 yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), yt(:knon, :), &
361                 u1lay(:knon), v1lay(:knon), coefh(:knon, :), yt, yq, &                 yq(:knon, :), yts(:knon), ypaprs(:knon, :), ypplay(:knon, :), &
362                 yts(:knon), ypaprs, ypplay, ydelp, yrads(:knon), yalb(:knon), &                 ydelp(:knon, :), radsol(:knon), yalb(:knon), snow(:knon), &
363                 snow(:knon), yqsurf, yrain_f, ysnow_f, yfluxlat(:knon), &                 yqsurf(:knon), yrain_fall(:knon), ysnow_fall(:knon), &
364                 pctsrf_new_sic, yagesno(:knon), y_d_t, y_d_q, y_d_ts(:knon), &                 yfluxlat(:knon), pctsrf_new_sic(ni(:knon)), yagesno(:knon), &
365                 yz0_new, y_flux_t(:knon), y_flux_q(:knon), y_dflux_t(:knon), &                 y_d_t(:knon, :), y_d_q(:knon, :), y_d_ts(:knon), &
366                 y_dflux_q(:knon), y_fqcalving, y_ffonte, y_run_off_lic_0)                 yz0_new(:knon), y_flux_t(:knon), y_flux_q(:knon), &
367                   y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving(:knon), &
368                   y_ffonte(:knon), y_run_off_lic_0(:knon), y_run_off_lic(:knon))
369    
370            ! calculer la longueur de rugosite sur ocean            ! calculer la longueur de rugosite sur ocean
371    
372            yrugm = 0.            yrugm = 0.
373    
374            IF (nsrf == is_oce) THEN            IF (nsrf == is_oce) THEN
375               DO j = 1, knon               DO j = 1, knon
376                  yrugm(j) = 0.018 * coefm(j, 1) * (u1lay(j)**2 + v1lay(j)**2) &                  yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
377                       / rg + 0.11 * 14E-6 &                       / rg + 0.11 * 14E-6 &
378                       / sqrt(coefm(j, 1) * (u1lay(j)**2 + v1lay(j)**2))                       / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
379                  yrugm(j) = max(1.5E-05, yrugm(j))                  yrugm(j) = max(1.5E-05, yrugm(j))
380               END DO               END DO
381            END IF            END IF
           DO j = 1, knon  
              y_dflux_t(j) = y_dflux_t(j) * ypct(j)  
              y_dflux_q(j) = y_dflux_q(j) * ypct(j)  
           END DO  
382    
383            DO k = 1, klev            DO k = 1, klev
384               DO j = 1, knon               DO j = 1, knon
385                  i = ni(j)                  i = ni(j)
386                  coefh(j, k) = coefh(j, k) * ypct(j)                  y_d_t(j, k) = y_d_t(j, k) * ypctsrf(j)
387                  coefm(j, k) = coefm(j, k) * ypct(j)                  y_d_q(j, k) = y_d_q(j, k) * ypctsrf(j)
388                  y_d_t(j, k) = y_d_t(j, k) * ypct(j)                  y_d_u(j, k) = y_d_u(j, k) * ypctsrf(j)
389                  y_d_q(j, k) = y_d_q(j, k) * ypct(j)                  y_d_v(j, k) = y_d_v(j, k) * ypctsrf(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)  
390               END DO               END DO
391            END DO            END DO
392    
# Line 450  contains Line 395  contains
395            flux_u(ni(:knon), nsrf) = y_flux_u(:knon)            flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
396            flux_v(ni(:knon), nsrf) = y_flux_v(:knon)            flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
397    
           evap(:, nsrf) = -flux_q(:, nsrf)  
   
398            falbe(:, nsrf) = 0.            falbe(:, nsrf) = 0.
399            fsnow(:, nsrf) = 0.            fsnow(:, nsrf) = 0.
400            qsurf(:, nsrf) = 0.            fqsurf(:, nsrf) = 0.
401            frugs(:, nsrf) = 0.            frugs(:, nsrf) = 0.
402            DO j = 1, knon            DO j = 1, knon
403               i = ni(j)               i = ni(j)
404               d_ts(i, nsrf) = y_d_ts(j)               d_ts(i, nsrf) = y_d_ts(j)
405               falbe(i, nsrf) = yalb(j)               falbe(i, nsrf) = yalb(j)
406               fsnow(i, nsrf) = snow(j)               fsnow(i, nsrf) = snow(j)
407               qsurf(i, nsrf) = yqsurf(j)               fqsurf(i, nsrf) = yqsurf(j)
408               frugs(i, nsrf) = yz0_new(j)               frugs(i, nsrf) = yz0_new(j)
409               fluxlat(i, nsrf) = yfluxlat(j)               fluxlat(i, nsrf) = yfluxlat(j)
410               IF (nsrf == is_oce) THEN               IF (nsrf == is_oce) THEN
# Line 471  contains Line 414  contains
414               agesno(i, nsrf) = yagesno(j)               agesno(i, nsrf) = yagesno(j)
415               fqcalving(i, nsrf) = y_fqcalving(j)               fqcalving(i, nsrf) = y_fqcalving(j)
416               ffonte(i, nsrf) = y_ffonte(j)               ffonte(i, nsrf) = y_ffonte(j)
417               cdragh(i) = cdragh(i) + coefh(j, 1)               cdragh(i) = cdragh(i) + ycdragh(j) * ypctsrf(j)
418               cdragm(i) = cdragm(i) + coefm(j, 1)               cdragm(i) = cdragm(i) + ycdragm(j) * ypctsrf(j)
419               dflux_t(i) = dflux_t(i) + y_dflux_t(j)               dflux_t(i) = dflux_t(i) + y_dflux_t(j) * ypctsrf(j)
420               dflux_q(i) = dflux_q(i) + y_dflux_q(j)               dflux_q(i) = dflux_q(i) + y_dflux_q(j) * ypctsrf(j)
              zu1(i) = zu1(i) + u1lay(j) * ypct(j)  
              zv1(i) = zv1(i) + v1lay(j) * ypct(j)  
421            END DO            END DO
422            IF (nsrf == is_ter) THEN            IF (nsrf == is_ter) THEN
423               qsol(ni(:knon)) = yqsol(:knon)               qsol(ni(:knon)) = yqsol(:knon)
# Line 484  contains Line 425  contains
425               DO j = 1, knon               DO j = 1, knon
426                  i = ni(j)                  i = ni(j)
427                  run_off_lic_0(i) = y_run_off_lic_0(j)                  run_off_lic_0(i) = y_run_off_lic_0(j)
428                    run_off_lic(i) = y_run_off_lic(j)
429               END DO               END DO
430            END IF            END IF
431    
# Line 497  contains Line 439  contains
439                  d_q(i, k) = d_q(i, k) + y_d_q(j, k)                  d_q(i, k) = d_q(i, k) + y_d_q(j, k)
440                  d_u(i, k) = d_u(i, k) + y_d_u(j, k)                  d_u(i, k) = d_u(i, k) + y_d_u(j, k)
441                  d_v(i, k) = d_v(i, k) + y_d_v(j, k)                  d_v(i, k) = d_v(i, k) + y_d_v(j, k)
                 ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)  
442               END DO               END DO
443            END DO            END DO
444    
445              forall (k = 2:klev) coefh(ni(:knon), k) &
446                   = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypctsrf(:knon)
447    
448            ! diagnostic t, q a 2m et u, v a 10m            ! diagnostic t, q a 2m et u, v a 10m
449    
450            DO j = 1, knon            DO j = 1, knon
451               i = ni(j)               i = ni(j)
452               uzon(j) = yu(j, 1) + y_d_u(j, 1)               u1(j) = yu(j, 1) + y_d_u(j, 1)
453               vmer(j) = yv(j, 1) + y_d_v(j, 1)               v1(j) = yv(j, 1) + y_d_v(j, 1)
454               tair1(j) = yt(j, 1) + y_d_t(j, 1)               tair1(j) = yt(j, 1) + y_d_t(j, 1)
455               qair1(j) = yq(j, 1) + y_d_q(j, 1)               qair1(j) = yq(j, 1) + y_d_q(j, 1)
456               zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &               zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
# Line 518  contains Line 462  contains
462               END IF               END IF
463               psfce(j) = ypaprs(j, 1)               psfce(j) = ypaprs(j, 1)
464               patm(j) = ypplay(j, 1)               patm(j) = ypplay(j, 1)
   
              qairsol(j) = yqsurf(j)  
465            END DO            END DO
466    
467            CALL stdlevvar(klon, knon, nsrf, zxli, uzon(:knon), vmer(:knon), &            CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), qair1, &
468                 tair1, qair1, zgeo1, tairsol, qairsol, rugo1, psfce, patm, &                 zgeo1, tairsol, yqsurf(:knon), rugo1, psfce, patm, yt2m, yq2m, &
469                 yt2m, yq2m, yt10m, yq10m, yu10m, yustar)                 yt10m, yq10m, wind10m(:knon), ustar(:knon))
470    
471            DO j = 1, knon            DO j = 1, knon
472               i = ni(j)               i = ni(j)
473               t2m(i, nsrf) = yt2m(j)               t2m(i, nsrf) = yt2m(j)
474               q2m(i, nsrf) = yq2m(j)               q2m(i, nsrf) = yq2m(j)
475    
476               u10m_srf(i, nsrf) = (yu10m(j) * uzon(j)) &               u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
477                    / sqrt(uzon(j)**2 + vmer(j)**2)                    / sqrt(u1(j)**2 + v1(j)**2)
478               v10m_srf(i, nsrf) = (yu10m(j) * vmer(j)) &               v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
479                    / sqrt(uzon(j)**2 + vmer(j)**2)                    / sqrt(u1(j)**2 + v1(j)**2)
480            END DO            END DO
481    
482            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, yustar, y_flux_t(:knon), &            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
483                 y_flux_q(:knon), yu, yv, yt, yq, ypblh(:knon), ycapcl, &                 y_flux_q(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
484                 yoliqcl, ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)                 yq(:knon, :), ypblh(:knon), ycapcl, yoliqcl, ycteicl, ypblt, &
485                   ytherm, ylcl)
486    
487            DO j = 1, knon            DO j = 1, knon
488               i = ni(j)               i = ni(j)
# Line 550  contains Line 493  contains
493               cteicl(i, nsrf) = ycteicl(j)               cteicl(i, nsrf) = ycteicl(j)
494               pblt(i, nsrf) = ypblt(j)               pblt(i, nsrf) = ypblt(j)
495               therm(i, nsrf) = ytherm(j)               therm(i, nsrf) = ytherm(j)
              trmb1(i, nsrf) = ytrmb1(j)  
              trmb2(i, nsrf) = ytrmb2(j)  
              trmb3(i, nsrf) = ytrmb3(j)  
496            END DO            END DO
497    
498            DO j = 1, knon            IF (iflag_pbl >= 6) q2(ni(:knon), :, nsrf) = yq2(:knon, :)
              DO k = 1, klev + 1  
                 i = ni(j)  
                 q2(i, k, nsrf) = yq2(j, k)  
              END DO  
           END DO  
499         else         else
500            fsnow(:, nsrf) = 0.            fsnow(:, nsrf) = 0.
501         end IF if_knon         end IF if_knon
# Line 571  contains Line 506  contains
506      pctsrf(:, is_oce) = pctsrf_new_oce      pctsrf(:, is_oce) = pctsrf_new_oce
507      pctsrf(:, is_sic) = pctsrf_new_sic      pctsrf(:, is_sic) = pctsrf_new_sic
508    
509      firstcal = .false.      CALL histwrite_phy("run_off_lic", run_off_lic)
510    
511    END SUBROUTINE clmain    END SUBROUTINE pbl_surface
512    
513  end module clmain_m  end module pbl_surface_m

Legend:
Removed from v.225  
changed lines
  Added in v.324

  ViewVC Help
Powered by ViewVC 1.1.21