/[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 239 by guez, Fri Nov 10 15:16:48 2017 UTC trunk/phylmd/Interface_surf/pbl_surface.f90 revision 343 by guez, Mon Oct 28 08:14:26 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_seri, q_seri, u_seri, v_seri, julien, mu0, &
8         cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, pplay, fsnow, &         ftsol, cdmmax, cdhmax, ftsoil, qsol, paprs, play, fsnow, fqsurf, falbe, &
9         qsurf, evap, falbe, fluxlat, rain_fall, snow_f, fsolsw, fsollw, frugs, &         fluxlat, rain_fall, snow_fall, frugs, agesno, rugoro, d_t, d_q, d_u, &
10         agesno, rugoro, d_t, d_q, d_u, d_v, d_ts, flux_t, flux_q, flux_u, &         d_v, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, dflux_t, &
11         flux_v, cdragh, cdragm, q2, dflux_t, dflux_q, ycoefh, t2m, q2m, &         dflux_q, coefh, t2m, q2m, u10m_srf, v10m_srf, pblh, capcl, oliqcl, &
12         u10m_srf, v10m_srf, pblh, capcl, oliqcl, cteicl, pblt, therm, trmb1, &         cteicl, pblt, therm, plcl, fqcalving, ffonte, run_off_lic_0, albsol, &
13         trmb2, trmb3, plcl, fqcalving, ffonte, run_off_lic_0)         sollw, solsw, 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        use cdrag_m, only: cdrag
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
     use coefkzmin_m, only: coefkzmin  
     use coefkz2_m, only: coefkz2  
29      USE conf_gcm_m, ONLY: lmt_pas      USE conf_gcm_m, ONLY: 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 yamada4_m, only: yamada4  
41    
42      REAL, INTENT(IN):: dtime ! interval du temps (secondes)      REAL, INTENT(inout):: pctsrf(:, :) ! (klon, nbsrf)
43        ! pourcentages de surface de chaque maille
44    
45      REAL, INTENT(inout):: pctsrf(klon, nbsrf)      REAL, INTENT(IN):: t_seri(:, :) ! (klon, klev) air temperature, in K
46      ! tableau des pourcentages de surface de chaque maille      REAL, INTENT(IN):: q_seri(:, :) ! (klon, klev) mass fraction of water vapor
47        REAL, INTENT(IN):: u_seri(:, :), v_seri(:, :) ! (klon, klev) wind, in m s -1
     REAL, INTENT(IN):: t(klon, klev) ! temperature (K)  
     REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg / kg)  
     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(INout):: 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 61  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  
65      REAL qsurf(klon, nbsrf)      REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf)
66      REAL evap(klon, nbsrf)      ! column-density of mass of snow at the surface, in kg m-2
67    
68        REAL, INTENT(inout):: fqsurf(klon, nbsrf)
69      REAL, intent(inout):: falbe(klon, nbsrf)      REAL, intent(inout):: falbe(klon, nbsrf)
70    
71      REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)      REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)
72        ! flux de chaleur latente, en W m-2
73    
74      REAL, intent(in):: rain_fall(klon)      REAL, intent(in):: rain_fall(klon)
75      ! liquid water mass flux (kg / m2 / s), positive down      ! liquid water mass flux (kg / m2 / s), positive down
76    
77      REAL, intent(in):: snow_f(klon)      REAL, intent(in):: snow_fall(klon)
78      ! solid water mass flux (kg / m2 / s), positive down      ! solid water mass flux (kg / m2 / s), positive down
79    
     REAL, INTENT(IN):: fsolsw(klon, nbsrf), fsollw(klon, nbsrf)  
80      REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)      REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)
81      real agesno(klon, nbsrf)      real, intent(inout):: agesno(:, :) ! (klon, nbsrf)
82      REAL, INTENT(IN):: rugoro(klon)      REAL, INTENT(IN):: rugoro(klon)
83    
84      REAL d_t(klon, klev), d_q(klon, klev)      REAL, intent(out):: d_t(:, :), d_q(:, :) ! (klon, klev)
85      ! d_t------output-R- le changement pour "t"      ! changement pour t_seri et q_seri
     ! d_q------output-R- le changement pour "q"  
86    
87      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
88      ! changement pour "u" et "v"      ! changement pour "u_seri" et "v_seri"
   
     REAL, intent(out):: d_ts(:, :) ! (klon, nbsrf) variation of ftsol  
89    
90      REAL, intent(out):: flux_t(klon, nbsrf)      REAL, intent(out):: flux_t(klon, nbsrf)
91      ! flux de chaleur sensible (Cp T) (W / m2) (orientation positive vers      ! flux de chaleur sensible (c_p T) (W / m2) (orientation positive
92      ! le bas) à la surface      ! vers le bas) à la surface
93    
94      REAL, intent(out):: flux_q(klon, nbsrf)      REAL, intent(out):: flux_q(klon, nbsrf)
95      ! flux de vapeur d'eau (kg / m2 / s) à la surface      ! flux de vapeur d'eau (kg / m2 / s) à la surface
96    
97      REAL, intent(out):: flux_u(klon, nbsrf), flux_v(klon, nbsrf)      REAL, intent(out):: flux_u(:, :), flux_v(:, :) ! (klon, nbsrf)
98      ! tension du vent (flux turbulent de vent) à la surface, en Pa      ! tension du vent (flux turbulent de vent) à la surface, en Pa
99    
100      REAL, INTENT(out):: cdragh(klon), cdragm(klon)      REAL, INTENT(out):: cdragh(klon), cdragm(klon)
101      real q2(klon, klev + 1, nbsrf)      real q2(klon, klev + 1, nbsrf)
102    
103      REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)      ! Ocean slab:
104      ! dflux_t derive du flux sensible      REAL, INTENT(out):: dflux_t(klon) ! derive du flux sensible
105      ! dflux_q derive du flux latent      REAL, INTENT(out):: dflux_q(klon) ! derive du flux latent
     ! IM "slab" ocean  
106    
107      REAL, intent(out):: ycoefh(:, :) ! (klon, klev)      REAL, intent(out):: coefh(:, 2:) ! (klon, 2:klev)
108      ! Pour pouvoir extraire les coefficients d'\'echange, le champ      ! Pour pouvoir extraire les coefficients d'\'echange, le champ
109      ! "ycoefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de      ! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
110      ! ce champ sur les quatre sous-surfaces du mod\`ele.      ! ce champ sur les quatre sous-surfaces du mod\`ele.
111    
112      REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)      REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)
# Line 123  contains Line 121  contains
121      REAL capcl(klon, nbsrf)      REAL capcl(klon, nbsrf)
122      REAL oliqcl(klon, nbsrf)      REAL oliqcl(klon, nbsrf)
123      REAL cteicl(klon, nbsrf)      REAL cteicl(klon, nbsrf)
124      REAL, INTENT(inout):: pblt(klon, nbsrf) ! T au nveau HCL      REAL, INTENT(inout):: pblt(:, :) ! (klon, nbsrf) temp\'erature au nveau HCL
125      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  
126      REAL plcl(klon, nbsrf)      REAL plcl(klon, nbsrf)
127      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)  
128      ! ffonte----Flux thermique utilise pour fondre la neige      REAL, intent(out):: fqcalving(klon, nbsrf)
129      ! 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
130      !           hauteur de neige, en kg / m2 / s      ! hauteur de neige, en kg / m2 / s
131      REAL run_off_lic_0(klon)  
132        real ffonte(klon, nbsrf) ! flux thermique utilise pour fondre la neige
133        REAL, intent(inout):: run_off_lic_0(:) ! (klon)
134    
135        REAL, intent(out):: albsol(:) ! (klon)
136        ! albedo du sol total, visible, moyen par maille
137    
138        REAL, intent(in):: sollw(:) ! (klon)
139        ! surface net downward longwave flux, in W m-2
140    
141        REAL, intent(in):: solsw(:) ! (klon)
142        ! surface net downward shortwave flux, in W m-2
143    
144        REAL, intent(in):: tsol(:) ! (klon)
145    
146      ! Local:      ! Local:
147    
148      LOGICAL:: firstcal = .true.      REAL d_ts(klon, nbsrf) ! variation of ftsol
149        REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
150        REAL fsolsw(klon, nbsrf) ! flux solaire absorb\'e pour chaque sous-surface
151    
152      ! la nouvelle repartition des surfaces sortie de l'interface      ! la nouvelle repartition des surfaces sortie de l'interface
153      REAL, save:: pctsrf_new_oce(klon)      REAL, save:: pctsrf_new_oce(klon)
154      REAL, save:: pctsrf_new_sic(klon)      REAL, save:: pctsrf_new_sic(klon)
155    
156      REAL y_fqcalving(klon), y_ffonte(klon)      REAL y_fqcalving(klon), y_ffonte(klon)
157      real y_run_off_lic_0(klon)      real y_run_off_lic_0(klon), y_run_off_lic(klon)
158        REAL run_off_lic(klon) ! ruissellement total
159      REAL rugmer(klon)      REAL rugmer(klon)
160      REAL ytsoil(klon, nsoilmx)      REAL ytsoil(klon, nsoilmx)
161      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      REAL yts(klon), ypctsrf(klon), yz0_new(klon)
162        real yrugos(klon) ! longueur de rugosit\'e, en m
163      REAL yalb(klon)      REAL yalb(klon)
164      REAL snow(klon), yqsurf(klon), yagesno(klon)      REAL snow(klon) ! column-density of mass of snow at the surface, in kg m-2
165        real yqsurf(klon), yagesno(klon)
166      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
167      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
168      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
169      REAL yrugm(klon), yrads(klon), yrugoro(klon)      REAL yrugm(klon), radsol(klon), yrugoro(klon)
170      REAL yfluxlat(klon)      REAL yfluxlat(klon)
171      REAL y_d_ts(klon)      REAL y_d_ts(klon)
172      REAL y_d_t(klon, klev), y_d_q(klon, klev)      REAL y_d_t(klon, klev), y_d_q(klon, klev)
# Line 164  contains Line 174  contains
174      REAL y_flux_t(klon), y_flux_q(klon)      REAL y_flux_t(klon), y_flux_q(klon)
175      REAL y_flux_u(klon), y_flux_v(klon)      REAL y_flux_u(klon), y_flux_v(klon)
176      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
177      REAL coefh(klon, 2:klev), coefm(klon, 2:klev)      REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
178      real ycdragh(klon), ycdragm(klon)      real ycdragh(klon), ycdragm(klon)
179      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
180      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
181      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), zlev(klon, klev + 1), yteta(klon, klev)  
     REAL ykmm(klon, klev + 1), ykmn(klon, klev + 1)  
182      REAL yq2(klon, klev + 1)      REAL yq2(klon, klev + 1)
183      REAL delp(klon, klev)      REAL delp(klon, klev)
184      INTEGER i, k, nsrf      INTEGER i, k, nsrf
# Line 192  contains Line 199  contains
199      REAL ycteicl(klon)      REAL ycteicl(klon)
200      REAL ypblt(klon)      REAL ypblt(klon)
201      REAL ytherm(klon)      REAL ytherm(klon)
     REAL ytrmb1(klon)  
     REAL ytrmb2(klon)  
     REAL ytrmb3(klon)  
202      REAL u1(klon), v1(klon)      REAL u1(klon), v1(klon)
203      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon)
     REAL psfce(klon), patm(klon)  
   
     REAL qairsol(klon), zgeo1(klon)  
204      REAL rugo1(klon)      REAL rugo1(klon)
205        REAL zgeop(klon, klev)
206    
207      !------------------------------------------------------------      !------------------------------------------------------------
208    
209        albsol = sum(falbe * pctsrf, dim = 2)
210    
211        ! R\'epartition sous maille des flux longwave et shortwave
212        ! R\'epartition du longwave par sous-surface lin\'earis\'ee
213    
214        forall (nsrf = 1:nbsrf)
215           fsollw(:, nsrf) = sollw + 4. * RSIGMA * tsol**3 &
216                * (tsol - ftsol(:, nsrf))
217           fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
218        END forall
219    
220      ytherm = 0.      ytherm = 0.
221    
222      DO k = 1, klev ! epaisseur de couche      DO k = 1, klev ! epaisseur de couche
# Line 218  contains Line 231  contains
231      cdragm = 0.      cdragm = 0.
232      dflux_t = 0.      dflux_t = 0.
233      dflux_q = 0.      dflux_q = 0.
     ypct = 0.  
     yqsurf = 0.  
     yrain_f = 0.  
     ysnow_f = 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 243  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 262  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 275  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               yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)            y_run_off_lic_0(:knon) = run_off_lic_0(ni(:knon))
              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 301  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_seri(i, k)
310                  yv(j, k) = v(i, k)                  yv(j, k) = v_seri(i, k)
311                  yt(j, k) = t(i, k)                  yt(j, k) = t_seri(i, k)
312                  yq(j, k) = q(i, k)                  yq(j, k) = q_seri(i, k)
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:
317            CALL coefkz(nsrf, ypaprs, ypplay, ksta, ksta_ter, yts(:knon), &  
318                 yrugos, yu, yv, yt, yq, yqsurf(:knon), coefm(:knon, :), &            zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
319                 coefh(:knon, :), ycdragm(:knon), ycdragh(:knon))                 + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
320    
321              DO k = 2, klev
322                 zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
323                      * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
324                      * (ypplay(:knon, k - 1) - ypplay(:knon, k))
325              ENDDO
326    
327              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 == 1) THEN            IF (iflag_pbl == 1) THEN
              CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0(:knon, 2:), &  
                   ycoefh0(:knon, 2:))  
              ycoefm0(:knon, 1) = 0.  
              ycoefh0(:knon, 1) = 0.  
              coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, 2:))  
              coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, 2:))  
333               ycdragm(:knon) = max(ycdragm(:knon), 0.)               ycdragm(:knon) = max(ycdragm(:knon), 0.)
334               ycdragh(:knon) = max(ycdragh(:knon), 0.)               ycdragh(:knon) = max(ycdragh(:knon), 0.)
335            END IF            end IF
336    
           ! on met un seuil pour ycdragm et ycdragh  
337            IF (nsrf == is_oce) THEN            IF (nsrf == is_oce) THEN
338                 ! On met un seuil pour ycdragm et ycdragh :
339               ycdragm(:knon) = min(ycdragm(:knon), cdmmax)               ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
340               ycdragh(:knon) = min(ycdragh(:knon), cdhmax)               ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
341            END IF            END IF
342    
343            IF (ok_kzmin) THEN            IF (iflag_pbl >= 6) yq2(:knon, :) = q2(ni(:knon), :, nsrf)
344               ! Calcul d'une diffusion minimale pour les conditions tres stables            call coef_diff_turb(nsrf, ni(:knon), ypaprs(:knon, :), &
345               CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &                 ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
346                    ycdragm(:knon), ycoefh0(:knon, 2:))                 yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
347               ycoefm0(:knon, 2:) = ycoefh0(:knon, 2:)                 ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
              coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, 2:))  
              coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, 2:))  
           END IF  
   
           IF (iflag_pbl >= 6) THEN  
              ! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et  
              ! Fr\'ed\'eric Hourdin  
              yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &  
                   + ypplay(:knon, 1))) &  
                   * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg  
   
              DO k = 2, klev  
                 yzlay(:knon, k) = yzlay(: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  
348    
349               zlev(:knon, 1) = 0.            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
              zlev(:knon, klev + 1) = 2. * yzlay(:knon, klev) &  
                   - yzlay(:knon, klev - 1)  
   
              DO k = 2, klev  
                 zlev(:knon, k) = 0.5 * (yzlay(:knon, k) + yzlay(: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  
   
              ustar(:knon) = ustarhb(yu(:knon, 1), yv(:knon, 1), ycdragm(:knon))  
              CALL yamada4(dtime, rg, zlev(:knon, :), yzlay(:knon, :), &  
                   yu(:knon, :), yv(:knon, :), yteta(:knon, :), yq2(:knon, :), &  
                   ykmm(:knon, :), ykmn(:knon, :), ustar(:knon))  
              coefm(:knon, :) = ykmm(:knon, 2:klev)  
              coefh(:knon, :) = ykmn(:knon, 2:klev)  
           END IF  
   
           CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), coefm(:knon, :), &  
350                 ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &                 ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
351                 ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &                 ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
352                 y_flux_u(:knon))                 y_flux_u(:knon))
353            CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), coefm(:knon, :), &            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
354                 ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &                 ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
355                 ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &                 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                 yu(:knon, 1), yv(:knon, 1), coefh(:knon, :), ycdragh(:knon), &                 yq(:knon, :), yts(:knon), ypaprs(:knon, :), ypplay(:knon, :), &
362                 yt, yq, yts(:knon), ypaprs, ypplay, ydelp, yrads(:knon), &                 ydelp(:knon, :), radsol(:knon), yalb(:knon), snow(:knon), &
363                 yalb(:knon), snow(:knon), yqsurf, yrain_f, ysnow_f, &                 yqsurf(:knon), yrain_fall(:knon), ysnow_fall(:knon), &
364                 yfluxlat(:knon), pctsrf_new_sic, yagesno(:knon), y_d_t, y_d_q, &                 yfluxlat(:knon), pctsrf_new_sic(ni(:knon)), yagesno(:knon), &
365                 y_d_ts(:knon), yz0_new, y_flux_t(:knon), y_flux_q(:knon), &                 y_d_t(:knon, :), y_d_q(:knon, :), y_d_ts(:knon), &
366                 y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving, y_ffonte, &                 yz0_new(:knon), y_flux_t(:knon), y_flux_q(:knon), &
367                 y_run_off_lic_0)                 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 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &                  yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
# Line 413  contains Line 379  contains
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    
           DO k = 2, klev  
              DO j = 1, knon  
                 i = ni(j)  
                 coefh(j, k) = coefh(j, k) * ypct(j)  
                 coefm(j, k) = coefm(j, k) * ypct(j)  
              END DO  
           END DO  
           DO j = 1, knon  
              i = ni(j)  
              ycdragh(j) = ycdragh(j) * ypct(j)  
              ycdragm(j) = ycdragm(j) * ypct(j)  
           END DO  
383            DO k = 1, klev            DO k = 1, klev
384               DO j = 1, knon               DO j = 1, knon
385                  i = ni(j)                  y_d_t(j, k) = y_d_t(j, k) * ypctsrf(j)
386                  y_d_t(j, k) = y_d_t(j, k) * ypct(j)                  y_d_q(j, k) = y_d_q(j, k) * ypctsrf(j)
387                  y_d_q(j, k) = y_d_q(j, k) * ypct(j)                  y_d_u(j, k) = y_d_u(j, k) * ypctsrf(j)
388                  y_d_u(j, k) = y_d_u(j, k) * ypct(j)                  y_d_v(j, k) = y_d_v(j, k) * ypctsrf(j)
                 y_d_v(j, k) = y_d_v(j, k) * ypct(j)  
389               END DO               END DO
390            END DO            END DO
391    
# Line 445  contains Line 394  contains
394            flux_u(ni(:knon), nsrf) = y_flux_u(:knon)            flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
395            flux_v(ni(:knon), nsrf) = y_flux_v(:knon)            flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
396    
           evap(:, nsrf) = -flux_q(:, nsrf)  
   
397            falbe(:, nsrf) = 0.            falbe(:, nsrf) = 0.
398            fsnow(:, nsrf) = 0.            fsnow(:, nsrf) = 0.
399            qsurf(:, nsrf) = 0.            fqsurf(:, nsrf) = 0.
400            frugs(:, nsrf) = 0.            frugs(:, nsrf) = 0.
401            DO j = 1, knon            DO j = 1, knon
402               i = ni(j)               i = ni(j)
403               d_ts(i, nsrf) = y_d_ts(j)               d_ts(i, nsrf) = y_d_ts(j)
404               falbe(i, nsrf) = yalb(j)               falbe(i, nsrf) = yalb(j)
405               fsnow(i, nsrf) = snow(j)               fsnow(i, nsrf) = snow(j)
406               qsurf(i, nsrf) = yqsurf(j)               fqsurf(i, nsrf) = yqsurf(j)
407               frugs(i, nsrf) = yz0_new(j)               frugs(i, nsrf) = yz0_new(j)
408               fluxlat(i, nsrf) = yfluxlat(j)               fluxlat(i, nsrf) = yfluxlat(j)
409               IF (nsrf == is_oce) THEN               IF (nsrf == is_oce) THEN
# Line 466  contains Line 413  contains
413               agesno(i, nsrf) = yagesno(j)               agesno(i, nsrf) = yagesno(j)
414               fqcalving(i, nsrf) = y_fqcalving(j)               fqcalving(i, nsrf) = y_fqcalving(j)
415               ffonte(i, nsrf) = y_ffonte(j)               ffonte(i, nsrf) = y_ffonte(j)
416               cdragh(i) = cdragh(i) + ycdragh(j)               cdragh(i) = cdragh(i) + ycdragh(j) * ypctsrf(j)
417               cdragm(i) = cdragm(i) + ycdragm(j)               cdragm(i) = cdragm(i) + ycdragm(j) * ypctsrf(j)
418               dflux_t(i) = dflux_t(i) + y_dflux_t(j)               dflux_t(i) = dflux_t(i) + y_dflux_t(j) * ypctsrf(j)
419               dflux_q(i) = dflux_q(i) + y_dflux_q(j)               dflux_q(i) = dflux_q(i) + y_dflux_q(j) * ypctsrf(j)
420            END DO            END DO
421            IF (nsrf == is_ter) THEN            IF (nsrf == is_ter) THEN
422               qsol(ni(:knon)) = yqsol(:knon)               qsol(ni(:knon)) = yqsol(:knon)
# Line 477  contains Line 424  contains
424               DO j = 1, knon               DO j = 1, knon
425                  i = ni(j)                  i = ni(j)
426                  run_off_lic_0(i) = y_run_off_lic_0(j)                  run_off_lic_0(i) = y_run_off_lic_0(j)
427                    run_off_lic(i) = y_run_off_lic(j)
428               END DO               END DO
429            END IF            END IF
430    
# Line 492  contains Line 440  contains
440                  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)
441               END DO               END DO
442            END DO            END DO
             
           DO j = 1, knon  
              i = ni(j)  
              DO k = 2, klev  
                 ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)  
              END DO  
           END DO  
443    
444            DO j = 1, knon            forall (k = 2:klev) coefh(ni(:knon), k) &
445               i = ni(j)                 = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypctsrf(:knon)
              ycoefh(i, 1) = ycoefh(i, 1) + ycdragh(j)  
           END DO  
446    
447            ! diagnostic t, q a 2m et u, v a 10m            ! Diagnostic temp\'erature, q \`a 2 m et u, v \`a 10 m:
448    
449            DO j = 1, knon            u1(:knon) = yu(:knon, 1) + y_d_u(:knon, 1)
450               i = ni(j)            v1(:knon) = yv(:knon, 1) + y_d_v(:knon, 1)
451               u1(j) = yu(j, 1) + y_d_u(j, 1)            tair1(:knon) = yt(:knon, 1) + y_d_t(:knon, 1)
              v1(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 == is_oce) THEN  
                 rugo1(j) = frugs(i, nsrf)  
              END IF  
              psfce(j) = ypaprs(j, 1)  
              patm(j) = ypplay(j, 1)  
452    
453               qairsol(j) = yqsurf(j)            IF (nsrf == is_oce) THEN
454            END DO               rugo1(:knon) = frugs(ni(:knon), is_oce)
455              else
456                 rugo1(:knon) = yrugos(:knon)
457              END IF
458    
459            CALL stdlevvar(klon, knon, nsrf, u1(:knon), v1(:knon), tair1(:knon), &            CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), &
460                 qair1, zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, &                 yq(:knon, 1) + y_d_q(:knon, 1), rd * tair1(:knon) &
461                 yq2m, yt10m, yq10m, wind10m(:knon), ustar)                 / (0.5 * (ypaprs(:knon, 1) + ypplay(:knon, 1))) &
462                   * (ypaprs(:knon, 1) - ypplay(:knon, 1)), &
463                   yts(:knon) + y_d_ts(:knon), yqsurf(:knon), rugo1, &
464                   ypaprs(:knon, 1), ypplay(:knon, 1), yt2m, yq2m, yt10m, yq10m, &
465                   wind10m(:knon), ustar(:knon))
466    
467            DO j = 1, knon            DO j = 1, knon
468               i = ni(j)               i = ni(j)
# Line 542  contains Line 476  contains
476            END DO            END DO
477    
478            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
479                 y_flux_q(:knon), yu, yv, yt, yq, ypblh(:knon), ycapcl, &                 y_flux_q(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
480                 yoliqcl, ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)                 yq(:knon, :), ypblh(:knon), ycapcl, yoliqcl, ycteicl, ypblt, &
481                   ytherm, ylcl)
482    
483            DO j = 1, knon            DO j = 1, knon
484               i = ni(j)               i = ni(j)
# Line 554  contains Line 489  contains
489               cteicl(i, nsrf) = ycteicl(j)               cteicl(i, nsrf) = ycteicl(j)
490               pblt(i, nsrf) = ypblt(j)               pblt(i, nsrf) = ypblt(j)
491               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)  
492            END DO            END DO
493    
494            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  
495         else         else
496            fsnow(:, nsrf) = 0.            fsnow(:, nsrf) = 0.
497         end IF if_knon         end IF if_knon
# Line 575  contains Line 502  contains
502      pctsrf(:, is_oce) = pctsrf_new_oce      pctsrf(:, is_oce) = pctsrf_new_oce
503      pctsrf(:, is_sic) = pctsrf_new_sic      pctsrf(:, is_sic) = pctsrf_new_sic
504    
505      firstcal = .false.      CALL histwrite_phy("run_off_lic", run_off_lic)
506        ftsol = ftsol + d_ts ! update surface temperature
507        CALL histwrite_phy("dtsvdfo", d_ts(:, is_oce))
508        CALL histwrite_phy("dtsvdft", d_ts(:, is_ter))
509        CALL histwrite_phy("dtsvdfg", d_ts(:, is_lic))
510        CALL histwrite_phy("dtsvdfi", d_ts(:, is_sic))
511    
512    END SUBROUTINE clmain    END SUBROUTINE pbl_surface
513    
514  end module clmain_m  end module pbl_surface_m

Legend:
Removed from v.239  
changed lines
  Added in v.343

  ViewVC Help
Powered by ViewVC 1.1.21