/[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/Sources/phylmd/clmain.f revision 206 by guez, Tue Aug 30 12:52:46 2016 UTC trunk/phylmd/pbl_surface.f revision 281 by guez, Fri Jul 20 16:28:36 2018 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, jour, rmu0, ts, cdmmax, &    SUBROUTINE pbl_surface(dtime, pctsrf, t, q, u, v, julien, mu0, ftsol, &
8         cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, pplay, snow, &         cdmmax, cdhmax, ftsoil, qsol, paprs, pplay, fsnow, qsurf, evap, falbe, &
9         qsurf, evap, falbe, fluxlat, rain_fall, snow_f, solsw, sollw, fder, &         fluxlat, rain_fall, snow_f, fsolsw, fsollw, frugs, agesno, rugoro, d_t, &
10         rlat, rugos, agesno, rugoro, d_t, d_q, d_u, d_v, d_ts, flux_t, flux_q, &         d_q, d_u, d_v, d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, &
11         flux_u, flux_v, cdragh, cdragm, q2, dflux_t, dflux_q, ycoefh, zu1, &         q2, dflux_t, dflux_q, coefh, t2m, q2m, u10m_srf, v10m_srf, pblh, capcl, &
12         zv1, t2m, q2m, u10m, v10m, pblh, capcl, oliqcl, cteicl, pblt, therm, &         oliqcl, cteicl, pblt, therm, plcl, fqcalving, ffonte, run_off_lic_0)
        trmb1, trmb2, trmb3, plcl, fqcalving, ffonte, run_off_lic_0)  
13    
14      ! 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
15      ! Author: Z. X. Li (LMD/CNRS), date: 1993/08/18      ! Author: Z. X. Li (LMD/CNRS), date: 1993 Aug. 18th
16      ! Objet : interface de couche limite (diffusion verticale)      ! Objet : interface de couche limite (diffusion verticale)
17    
18      ! 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 20  contains
20      ! ne tient pas compte de la diff\'erentiation des sous-fractions      ! ne tient pas compte de la diff\'erentiation des sous-fractions
21      ! de sol.      ! de sol.
22    
23      ! 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.  
   
24      use clqh_m, only: clqh      use clqh_m, only: clqh
25      use clvent_m, only: clvent      use clvent_m, only: clvent
26      use coefkz_m, only: coefkz      use coef_diff_turb_m, only: coef_diff_turb
27      use coefkzmin_m, only: coefkzmin      USE conf_gcm_m, ONLY: lmt_pas
     USE conf_gcm_m, ONLY: prt_level, lmt_pas  
28      USE conf_phys_m, ONLY: iflag_pbl      USE conf_phys_m, ONLY: iflag_pbl
29      USE dimphy, ONLY: klev, klon, zmasq      USE dimphy, ONLY: klev, klon
30      USE dimsoil, ONLY: nsoilmx      USE dimsoil, ONLY: nsoilmx
31      use hbtm_m, only: hbtm      use hbtm_m, only: hbtm
32      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
33      USE interfoce_lim_m, ONLY: interfoce_lim      USE interfoce_lim_m, ONLY: interfoce_lim
34        use phyetat0_m, only: zmasq
35      use stdlevvar_m, only: stdlevvar      use stdlevvar_m, only: stdlevvar
36      USE suphec_m, ONLY: rd, rg, rkappa      USE suphec_m, ONLY: rd, rg
37      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  
38    
39      REAL, INTENT(IN):: dtime ! interval du temps (secondes)      REAL, INTENT(IN):: dtime ! interval du temps (secondes)
40    
# Line 50  contains Line 42  contains
42      ! tableau des pourcentages de surface de chaque maille      ! tableau des pourcentages de surface de chaque maille
43    
44      REAL, INTENT(IN):: t(klon, klev) ! temperature (K)      REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
45      REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg/kg)      REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg / kg)
46      REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse      REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
47      INTEGER, INTENT(IN):: jour ! jour de l'annee en cours      INTEGER, INTENT(IN):: julien ! jour de l'annee en cours
48      REAL, intent(in):: rmu0(klon) ! cosinus de l'angle solaire zenithal          REAL, intent(in):: mu0(klon) ! cosinus de l'angle solaire zenithal    
49      REAL, INTENT(IN):: ts(klon, nbsrf) ! temperature du sol (en Kelvin)      REAL, INTENT(IN):: ftsol(:, :) ! (klon, nbsrf) temp\'erature du sol (en K)
50      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  
51    
52      REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)      REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
53      ! soil temperature of surface fraction      ! soil temperature of surface fraction
54    
55      REAL, INTENT(inout):: qsol(klon)      REAL, INTENT(inout):: qsol(:) ! (klon)
56      ! column-density of water in soil, in kg m-2      ! column-density of water in soil, in kg m-2
57    
58      REAL, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)      REAL, INTENT(IN):: paprs(klon, klev + 1) ! pression a intercouche (Pa)
59      REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)      REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
60      REAL, INTENT(inout):: snow(klon, nbsrf)      REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse
61      REAL qsurf(klon, nbsrf)      REAL qsurf(klon, nbsrf)
62      REAL evap(klon, nbsrf)      REAL evap(klon, nbsrf)
63      REAL, intent(inout):: falbe(klon, nbsrf)      REAL, intent(inout):: falbe(klon, nbsrf)
64        REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)
     REAL fluxlat(klon, nbsrf)  
65    
66      REAL, intent(in):: rain_fall(klon)      REAL, intent(in):: rain_fall(klon)
67      ! liquid water mass flux (kg/m2/s), positive down      ! liquid water mass flux (kg / m2 / s), positive down
68    
69      REAL, intent(in):: snow_f(klon)      REAL, intent(in):: snow_f(klon)
70      ! solid water mass flux (kg/m2/s), positive down      ! solid water mass flux (kg / m2 / s), positive down
   
     REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)  
     REAL, intent(in):: fder(klon)  
     REAL, INTENT(IN):: rlat(klon) ! latitude en degr\'es  
   
     REAL, intent(inout):: rugos(klon, nbsrf) ! longueur de rugosit\'e (en m)  
71    
72        REAL, INTENT(IN):: fsolsw(klon, nbsrf), fsollw(klon, nbsrf)
73        REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)
74      real agesno(klon, nbsrf)      real agesno(klon, nbsrf)
75      REAL, INTENT(IN):: rugoro(klon)      REAL, INTENT(IN):: rugoro(klon)
76    
77      REAL d_t(klon, klev), d_q(klon, klev)      REAL, intent(out):: d_t(klon, klev), d_q(klon, klev)
78      ! d_t------output-R- le changement pour "t"      ! changement pour t et q
     ! d_q------output-R- le changement pour "q"  
79    
80      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
81      ! changement pour "u" et "v"      ! changement pour "u" et "v"
82    
83      REAL, intent(out):: d_ts(klon, nbsrf) ! le changement pour "ts"      REAL, intent(out):: d_ts(:, :) ! (klon, nbsrf) variation of ftsol
84    
85      REAL, intent(out):: flux_t(klon, nbsrf)      REAL, intent(out):: flux_t(klon, nbsrf)
86      ! flux de chaleur sensible (Cp T) (W/m2) (orientation positive vers      ! flux de chaleur sensible (Cp T) (W / m2) (orientation positive vers
87      ! le bas) à la surface      ! le bas) à la surface
88    
89      REAL, intent(out):: flux_q(klon, nbsrf)      REAL, intent(out):: flux_q(klon, nbsrf)
90      ! flux de vapeur d'eau (kg/m2/s) à la surface      ! flux de vapeur d'eau (kg / m2 / s) à la surface
91    
92      REAL, intent(out):: flux_u(klon, nbsrf), flux_v(klon, nbsrf)      REAL, intent(out):: flux_u(klon, nbsrf), flux_v(klon, nbsrf)
93      ! tension du vent à la surface, en Pa      ! tension du vent (flux turbulent de vent) à la surface, en Pa
94    
95      REAL, INTENT(out):: cdragh(klon), cdragm(klon)      REAL, INTENT(out):: cdragh(klon), cdragm(klon)
96      real q2(klon, klev+1, nbsrf)      real q2(klon, klev + 1, nbsrf)
97    
98      REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)      REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)
99      ! dflux_t derive du flux sensible      ! dflux_t derive du flux sensible
100      ! dflux_q derive du flux latent      ! dflux_q derive du flux latent
101      ! IM "slab" ocean      ! IM "slab" ocean
102    
103      REAL, intent(out):: ycoefh(klon, klev)      REAL, intent(out):: coefh(:, 2:) ! (klon, 2:klev)
104      REAL, intent(out):: zu1(klon)      ! Pour pouvoir extraire les coefficients d'\'echange, le champ
105      REAL zv1(klon)      ! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
106      REAL t2m(klon, nbsrf), q2m(klon, nbsrf)      ! ce champ sur les quatre sous-surfaces du mod\`ele.
107      REAL u10m(klon, nbsrf), v10m(klon, nbsrf)  
108        REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)
109      ! Ionela Musat cf. Anne Mathieu : planetary boundary layer, hbtm  
110      ! (Comme les autres diagnostics on cumule dans physiq ce qui      REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)
111      ! permet de sortir les grandeurs par sous-surface)      ! composantes du vent \`a 10m sans spirale d'Ekman
112    
113        ! Ionela Musat. Cf. Anne Mathieu : planetary boundary layer, hbtm.
114        ! Comme les autres diagnostics on cumule dans physiq ce qui permet
115        ! de sortir les grandeurs par sous-surface.
116      REAL pblh(klon, nbsrf) ! height of planetary boundary layer      REAL pblh(klon, nbsrf) ! height of planetary boundary layer
117      REAL capcl(klon, nbsrf)      REAL capcl(klon, nbsrf)
118      REAL oliqcl(klon, nbsrf)      REAL oliqcl(klon, nbsrf)
119      REAL cteicl(klon, nbsrf)      REAL cteicl(klon, nbsrf)
120      REAL pblt(klon, nbsrf)      REAL, INTENT(inout):: pblt(klon, nbsrf) ! T au nveau HCL
     ! pblT------- T au nveau HCL  
121      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  
122      REAL plcl(klon, nbsrf)      REAL plcl(klon, nbsrf)
123      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)  
124        REAL, intent(out):: fqcalving(klon, nbsrf)
125        ! flux d'eau "perdue" par la surface et necessaire pour limiter la
126        ! hauteur de neige, en kg / m2 / s
127    
128        real ffonte(klon, nbsrf)
129      ! ffonte----Flux thermique utilise pour fondre la neige      ! ffonte----Flux thermique utilise pour fondre la neige
     ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la  
     !           hauteur de neige, en kg/m2/s  
130      REAL run_off_lic_0(klon)      REAL run_off_lic_0(klon)
131    
132      ! Local:      ! Local:
# Line 155  contains Line 139  contains
139    
140      REAL y_fqcalving(klon), y_ffonte(klon)      REAL y_fqcalving(klon), y_ffonte(klon)
141      real y_run_off_lic_0(klon)      real y_run_off_lic_0(klon)
   
142      REAL rugmer(klon)      REAL rugmer(klon)
   
143      REAL ytsoil(klon, nsoilmx)      REAL ytsoil(klon, nsoilmx)
144        REAL yts(klon), ypct(klon), yz0_new(klon)
145      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      real yrugos(klon) ! longeur de rugosite (en m)
146      REAL yalb(klon)      REAL yalb(klon)
147      REAL yu1(klon), yv1(klon)      REAL snow(klon), yqsurf(klon), yagesno(klon)
148      ! on rajoute en output yu1 et yv1 qui sont les vents dans      real yqsol(klon) ! column-density of water in soil, in kg m-2
149      ! la premiere couche      REAL yrain_f(klon) ! liquid water mass flux (kg / m2 / s), positive down
150      REAL ysnow(klon), yqsurf(klon), yagesno(klon)      REAL ysnow_f(klon) ! solid water mass flux (kg / m2 / s), positive down
   
     real yqsol(klon)  
     ! column-density of water in soil, in kg m-2  
   
     REAL yrain_f(klon)  
     ! liquid water mass flux (kg/m2/s), positive down  
   
     REAL ysnow_f(klon)  
     ! solid water mass flux (kg/m2/s), positive down  
   
     REAL yfder(klon)  
151      REAL yrugm(klon), yrads(klon), yrugoro(klon)      REAL yrugm(klon), yrads(klon), yrugoro(klon)
   
152      REAL yfluxlat(klon)      REAL yfluxlat(klon)
   
153      REAL y_d_ts(klon)      REAL y_d_ts(klon)
154      REAL y_d_t(klon, klev), y_d_q(klon, klev)      REAL y_d_t(klon, klev), y_d_q(klon, klev)
155      REAL y_d_u(klon, klev), y_d_v(klon, klev)      REAL y_d_u(klon, klev), y_d_v(klon, klev)
156      REAL y_flux_t(klon), y_flux_q(klon)      REAL y_flux_t(klon), y_flux_q(klon)
157      REAL y_flux_u(klon), y_flux_v(klon)      REAL y_flux_u(klon), y_flux_v(klon)
158      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
159      REAL coefh(klon, klev), coefm(klon, klev)      REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
160        real ycdragh(klon), ycdragm(klon)
161      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
162      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
163      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)      REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
164        REAL yq2(klon, klev + 1)
     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)  
     REAL yq2(klon, klev+1)  
     REAL q2diag(klon, klev+1)  
   
     REAL u1lay(klon), v1lay(klon)  
165      REAL delp(klon, klev)      REAL delp(klon, klev)
166      INTEGER i, k, nsrf      INTEGER i, k, nsrf
   
167      INTEGER ni(klon), knon, j      INTEGER ni(klon), knon, j
168    
169      REAL pctsrf_pot(klon, nbsrf)      REAL pctsrf_pot(klon, nbsrf)
170      ! "pourcentage potentiel" pour tenir compte des \'eventuelles      ! "pourcentage potentiel" pour tenir compte des \'eventuelles
171      ! apparitions ou disparitions de la glace de mer      ! apparitions ou disparitions de la glace de mer
172    
173      REAL zx_alf1, zx_alf2 ! valeur ambiante par extrapolation      REAL yt2m(klon), yq2m(klon), wind10m(klon)
174        REAL ustar(klon)
     REAL yt2m(klon), yq2m(klon), yu10m(klon)  
     REAL yustar(klon)  
175    
176      REAL yt10m(klon), yq10m(klon)      REAL yt10m(klon), yq10m(klon)
177      REAL ypblh(klon)      REAL ypblh(klon)
# Line 223  contains Line 181  contains
181      REAL ycteicl(klon)      REAL ycteicl(klon)
182      REAL ypblt(klon)      REAL ypblt(klon)
183      REAL ytherm(klon)      REAL ytherm(klon)
184      REAL ytrmb1(klon)      REAL u1(klon), v1(klon)
     REAL ytrmb2(klon)  
     REAL ytrmb3(klon)  
     REAL uzon(klon), vmer(klon)  
185      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
186      REAL psfce(klon), patm(klon)      REAL psfce(klon), patm(klon)
187    
188      REAL qairsol(klon), zgeo1(klon)      REAL qairsol(klon), zgeo1(klon)
189      REAL rugo1(klon)      REAL rugo1(klon)
190        REAL zgeop(klon, klev)
     ! utiliser un jeu de fonctions simples                
     LOGICAL zxli  
     PARAMETER (zxli=.FALSE.)  
191    
192      !------------------------------------------------------------      !------------------------------------------------------------
193    
# Line 243  contains Line 195  contains
195    
196      DO k = 1, klev ! epaisseur de couche      DO k = 1, klev ! epaisseur de couche
197         DO i = 1, klon         DO i = 1, klon
198            delp(i, k) = paprs(i, k) - paprs(i, k+1)            delp(i, k) = paprs(i, k) - paprs(i, k + 1)
199         END DO         END DO
200      END DO      END DO
     DO i = 1, klon ! vent de la premiere couche  
        zx_alf1 = 1.0  
        zx_alf2 = 1.0 - zx_alf1  
        u1lay(i) = u(i, 1)*zx_alf1 + u(i, 2)*zx_alf2  
        v1lay(i) = v(i, 1)*zx_alf1 + v(i, 2)*zx_alf2  
     END DO  
201    
202      ! Initialization:      ! Initialization:
203      rugmer = 0.      rugmer = 0.
# Line 259  contains Line 205  contains
205      cdragm = 0.      cdragm = 0.
206      dflux_t = 0.      dflux_t = 0.
207      dflux_q = 0.      dflux_q = 0.
     zu1 = 0.  
     zv1 = 0.  
208      ypct = 0.      ypct = 0.
     yts = 0.  
     ysnow = 0.  
209      yqsurf = 0.      yqsurf = 0.
210      yrain_f = 0.      yrain_f = 0.
211      ysnow_f = 0.      ysnow_f = 0.
     yfder = 0.  
212      yrugos = 0.      yrugos = 0.
     yu1 = 0.  
     yv1 = 0.  
     yrads = 0.  
213      ypaprs = 0.      ypaprs = 0.
214      ypplay = 0.      ypplay = 0.
215      ydelp = 0.      ydelp = 0.
# Line 281  contains Line 219  contains
219      yq = 0.      yq = 0.
220      y_dflux_t = 0.      y_dflux_t = 0.
221      y_dflux_q = 0.      y_dflux_q = 0.
     ytsoil = 999999.  
222      yrugoro = 0.      yrugoro = 0.
223      d_ts = 0.      d_ts = 0.
     yfluxlat = 0.  
224      flux_t = 0.      flux_t = 0.
225      flux_q = 0.      flux_q = 0.
226      flux_u = 0.      flux_u = 0.
227      flux_v = 0.      flux_v = 0.
228        fluxlat = 0.
229      d_t = 0.      d_t = 0.
230      d_q = 0.      d_q = 0.
231      d_u = 0.      d_u = 0.
232      d_v = 0.      d_v = 0.
233      ycoefh = 0.      coefh = 0.
234        fqcalving = 0.
235    
236      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
237      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
# Line 306  contains Line 244  contains
244    
245      ! Tester si c'est le moment de lire le fichier:      ! Tester si c'est le moment de lire le fichier:
246      if (mod(itap - 1, lmt_pas) == 0) then      if (mod(itap - 1, lmt_pas) == 0) then
247         CALL interfoce_lim(jour, pctsrf_new_oce, pctsrf_new_sic)         CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
248      endif      endif
249    
250      ! Boucler sur toutes les sous-fractions du sol:      ! Boucler sur toutes les sous-fractions du sol:
# Line 328  contains Line 266  contains
266            DO j = 1, knon            DO j = 1, knon
267               i = ni(j)               i = ni(j)
268               ypct(j) = pctsrf(i, nsrf)               ypct(j) = pctsrf(i, nsrf)
269               yts(j) = ts(i, nsrf)               yts(j) = ftsol(i, nsrf)
270               ysnow(j) = snow(i, nsrf)               snow(j) = fsnow(i, nsrf)
271               yqsurf(j) = qsurf(i, nsrf)               yqsurf(j) = qsurf(i, nsrf)
272               yalb(j) = falbe(i, nsrf)               yalb(j) = falbe(i, nsrf)
273               yrain_f(j) = rain_fall(i)               yrain_f(j) = rain_fall(i)
274               ysnow_f(j) = snow_f(i)               ysnow_f(j) = snow_f(i)
275               yagesno(j) = agesno(i, nsrf)               yagesno(j) = agesno(i, nsrf)
276               yfder(j) = fder(i)               yrugos(j) = frugs(i, nsrf)
              yrugos(j) = rugos(i, nsrf)  
277               yrugoro(j) = rugoro(i)               yrugoro(j) = rugoro(i)
278               yu1(j) = u1lay(i)               yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)
279               yv1(j) = v1lay(i)               ypaprs(j, klev + 1) = paprs(i, klev + 1)
              yrads(j) = solsw(i, nsrf) + sollw(i, nsrf)  
              ypaprs(j, klev+1) = paprs(i, klev+1)  
280               y_run_off_lic_0(j) = run_off_lic_0(i)               y_run_off_lic_0(j) = run_off_lic_0(i)
281            END DO            END DO
282    
283            ! For continent, copy soil water content            ! For continent, copy soil water content
284            IF (nsrf == is_ter) THEN            IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
              yqsol(:knon) = qsol(ni(:knon))  
           ELSE  
              yqsol = 0.  
           END IF  
285    
286            DO k = 1, nsoilmx            ytsoil(:knon, :) = ftsoil(ni(:knon), :, nsrf)
              DO j = 1, knon  
                 i = ni(j)  
                 ytsoil(j, k) = ftsoil(i, k, nsrf)  
              END DO  
           END DO  
287    
288            DO k = 1, klev            DO k = 1, klev
289               DO j = 1, knon               DO j = 1, knon
# Line 372  contains Line 298  contains
298               END DO               END DO
299            END DO            END DO
300    
301            ! calculer Cdrag et les coefficients d'echange            ! Calculer les géopotentiels de chaque couche:
302            CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &  
303                 yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))            zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
304                   + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
305    
306              DO k = 2, klev
307                 zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
308                      * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
309                      * (ypplay(:knon, k - 1) - ypplay(:knon, k))
310              ENDDO
311    
312              CALL cdrag(nsrf, sqrt(yu(:knon, 1)**2 + yv(:knon, 1)**2), &
313                   yt(:knon, 1), yq(:knon, 1), zgeop(:knon, 1), ypaprs(:knon, 1), &
314                   yts(:knon), yqsurf(:knon), yrugos(:knon), ycdragm(:knon), &
315                   ycdragh(:knon))
316    
317            IF (iflag_pbl == 1) THEN            IF (iflag_pbl == 1) THEN
318               CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)               ycdragm(:knon) = max(ycdragm(:knon), 0.)
319               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))               ycdragh(:knon) = max(ycdragh(:knon), 0.)
320               coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))            end IF
           END IF  
321    
322            ! on met un seuil pour coefm et coefh            ! on met un seuil pour ycdragm et ycdragh
323            IF (nsrf == is_oce) THEN            IF (nsrf == is_oce) THEN
324               coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)               ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
325               coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)               ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
           END IF  
   
           IF (ok_kzmin) THEN  
              ! Calcul d'une diffusion minimale pour les conditions tres stables  
              CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &  
                   coefm(:knon, 1), ycoefm0, ycoefh0)  
              coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))  
              coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))  
326            END IF            END IF
327    
328            IF (iflag_pbl >= 3) THEN            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(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  
329               DO k = 1, klev + 1               DO k = 1, klev + 1
330                  DO j = 1, knon                  DO j = 1, knon
331                     i = ni(j)                     i = ni(j)
332                     yq2(j, k) = q2(i, k, nsrf)                     yq2(j, k) = q2(i, k, nsrf)
333                  END DO                  END DO
334               END DO               END DO
335              end IF
336    
337               CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)            call coef_diff_turb(dtime, nsrf, ni(:knon), ypaprs(:knon, :), &
338               IF (prt_level > 9) PRINT *, 'USTAR = ', yustar                 ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
339                   yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
340               ! iflag_pbl peut \^etre utilis\'e comme longueur de m\'elange                 ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
341    
342               IF (iflag_pbl >= 11) THEN            CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
343                  CALL vdif_kcay(knon, dtime, rg, ypaprs, yzlev, yzlay, yu, yv, &                 ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
344                       yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, yustar, &                 ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
345                       iflag_pbl)                 y_flux_u(:knon))
346               ELSE            CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
347                  CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &                 ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
348                       coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)                 ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
349               END IF                 y_flux_v(:knon))
   
              coefm(:knon, 2:) = ykmm(:knon, 2:klev)  
              coefh(:knon, 2:) = ykmn(:knon, 2:klev)  
           END IF  
   
           ! calculer la diffusion des vitesses "u" et "v"  
           CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &  
                ypplay, ydelp, y_d_u, y_flux_u(:knon))  
           CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &  
                ypplay, ydelp, y_d_v, y_flux_v(:knon))  
350    
351            ! calculer la diffusion de "q" et de "h"            ! calculer la diffusion de "q" et de "h"
352            CALL clqh(dtime, jour, firstcal, rlat, nsrf, ni(:knon), ytsoil, &            CALL clqh(dtime, julien, firstcal, nsrf, ni(:knon), &
353                 yqsol, rmu0, yrugos, yrugoro, yu1, yv1, coefh(:knon, :), yt, &                 ytsoil(:knon, :), yqsol(:knon), mu0, yrugos, yrugoro, &
354                 yq, yts, ypaprs, ypplay, ydelp, yrads, yalb(:knon), ysnow, &                 yu(:knon, 1), yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), &
355                 yqsurf, yrain_f, ysnow_f, yfder, yfluxlat, pctsrf_new_sic, &                 yt(:knon, :), yq(:knon, :), yts(:knon), ypaprs(:knon, :), &
356                 yagesno(:knon), y_d_t, y_d_q, y_d_ts(:knon), yz0_new, &                 ypplay(:knon, :), ydelp, yrads(:knon), yalb(:knon), &
357                 y_flux_t(:knon), y_flux_q(:knon), y_dflux_t, y_dflux_q, &                 snow(:knon), yqsurf, yrain_f, ysnow_f, yfluxlat(:knon), &
358                 y_fqcalving, y_ffonte, y_run_off_lic_0)                 pctsrf_new_sic, yagesno(:knon), y_d_t(:knon, :), &
359                   y_d_q(:knon, :), y_d_ts(:knon), yz0_new, y_flux_t(:knon), &
360                   y_flux_q(:knon), y_dflux_t(:knon), y_dflux_q(:knon), &
361                   y_fqcalving(:knon), y_ffonte, y_run_off_lic_0)
362    
363            ! calculer la longueur de rugosite sur ocean            ! calculer la longueur de rugosite sur ocean
364            yrugm = 0.            yrugm = 0.
365            IF (nsrf == is_oce) THEN            IF (nsrf == is_oce) THEN
366               DO j = 1, knon               DO j = 1, knon
367                  yrugm(j) = 0.018*coefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &                  yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
368                       0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))                       / rg + 0.11 * 14E-6 &
369                         / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
370                  yrugm(j) = max(1.5E-05, yrugm(j))                  yrugm(j) = max(1.5E-05, yrugm(j))
371               END DO               END DO
372            END IF            END IF
373            DO j = 1, knon            DO j = 1, knon
374               y_dflux_t(j) = y_dflux_t(j)*ypct(j)               y_dflux_t(j) = y_dflux_t(j) * ypct(j)
375               y_dflux_q(j) = y_dflux_q(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)  
376            END DO            END DO
377    
378            DO k = 1, klev            DO k = 1, klev
379               DO j = 1, knon               DO j = 1, knon
380                  i = ni(j)                  i = ni(j)
381                  coefh(j, k) = coefh(j, k)*ypct(j)                  y_d_t(j, k) = y_d_t(j, k) * ypct(j)
382                  coefm(j, k) = coefm(j, k)*ypct(j)                  y_d_q(j, k) = y_d_q(j, k) * ypct(j)
383                  y_d_t(j, k) = y_d_t(j, k)*ypct(j)                  y_d_u(j, k) = y_d_u(j, k) * ypct(j)
384                  y_d_q(j, k) = y_d_q(j, k)*ypct(j)                  y_d_v(j, k) = y_d_v(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)  
385               END DO               END DO
386            END DO            END DO
387    
388            DO j = 1, knon            flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
389               i = ni(j)            flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
390               flux_t(i, nsrf) = y_flux_t(j)            flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
391               flux_q(i, nsrf) = y_flux_q(j)            flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
              flux_u(i, nsrf) = y_flux_u(j)  
              flux_v(i, nsrf) = y_flux_v(j)  
           END DO  
392    
393            evap(:, nsrf) = -flux_q(:, nsrf)            evap(:, nsrf) = -flux_q(:, nsrf)
394    
395            falbe(:, nsrf) = 0.            falbe(:, nsrf) = 0.
396            snow(:, nsrf) = 0.            fsnow(:, nsrf) = 0.
397            qsurf(:, nsrf) = 0.            qsurf(:, nsrf) = 0.
398            rugos(:, nsrf) = 0.            frugs(:, nsrf) = 0.
           fluxlat(:, nsrf) = 0.  
399            DO j = 1, knon            DO j = 1, knon
400               i = ni(j)               i = ni(j)
401               d_ts(i, nsrf) = y_d_ts(j)               d_ts(i, nsrf) = y_d_ts(j)
402               falbe(i, nsrf) = yalb(j)               falbe(i, nsrf) = yalb(j)
403               snow(i, nsrf) = ysnow(j)               fsnow(i, nsrf) = snow(j)
404               qsurf(i, nsrf) = yqsurf(j)               qsurf(i, nsrf) = yqsurf(j)
405               rugos(i, nsrf) = yz0_new(j)               frugs(i, nsrf) = yz0_new(j)
406               fluxlat(i, nsrf) = yfluxlat(j)               fluxlat(i, nsrf) = yfluxlat(j)
407               IF (nsrf == is_oce) THEN               IF (nsrf == is_oce) THEN
408                  rugmer(i) = yrugm(j)                  rugmer(i) = yrugm(j)
409                  rugos(i, nsrf) = yrugm(j)                  frugs(i, nsrf) = yrugm(j)
410               END IF               END IF
411               agesno(i, nsrf) = yagesno(j)               agesno(i, nsrf) = yagesno(j)
412               fqcalving(i, nsrf) = y_fqcalving(j)               fqcalving(i, nsrf) = y_fqcalving(j)
413               ffonte(i, nsrf) = y_ffonte(j)               ffonte(i, nsrf) = y_ffonte(j)
414               cdragh(i) = cdragh(i) + coefh(j, 1)               cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
415               cdragm(i) = cdragm(i) + coefm(j, 1)               cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
416               dflux_t(i) = dflux_t(i) + y_dflux_t(j)               dflux_t(i) = dflux_t(i) + y_dflux_t(j)
417               dflux_q(i) = dflux_q(i) + y_dflux_q(j)               dflux_q(i) = dflux_q(i) + y_dflux_q(j)
              zu1(i) = zu1(i) + yu1(j)  
              zv1(i) = zv1(i) + yv1(j)  
418            END DO            END DO
419            IF (nsrf == is_ter) THEN            IF (nsrf == is_ter) THEN
420               qsol(ni(:knon)) = yqsol(:knon)               qsol(ni(:knon)) = yqsol(:knon)
# Line 532  contains Line 426  contains
426            END IF            END IF
427    
428            ftsoil(:, :, nsrf) = 0.            ftsoil(:, :, nsrf) = 0.
429            DO k = 1, nsoilmx            ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)
              DO j = 1, knon  
                 i = ni(j)  
                 ftsoil(i, k, nsrf) = ytsoil(j, k)  
              END DO  
           END DO  
430    
431            DO j = 1, knon            DO j = 1, knon
432               i = ni(j)               i = ni(j)
# Line 546  contains Line 435  contains
435                  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)
436                  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)
437                  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)  
438               END DO               END DO
439            END DO            END DO
440    
441              forall (k = 2:klev) coefh(ni(:knon), k) &
442                   = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypct(:knon)
443    
444            ! diagnostic t, q a 2m et u, v a 10m            ! diagnostic t, q a 2m et u, v a 10m
445    
446            DO j = 1, knon            DO j = 1, knon
447               i = ni(j)               i = ni(j)
448               uzon(j) = yu(j, 1) + y_d_u(j, 1)               u1(j) = yu(j, 1) + y_d_u(j, 1)
449               vmer(j) = yv(j, 1) + y_d_v(j, 1)               v1(j) = yv(j, 1) + y_d_v(j, 1)
450               tair1(j) = yt(j, 1) + y_d_t(j, 1)               tair1(j) = yt(j, 1) + y_d_t(j, 1)
451               qair1(j) = yq(j, 1) + y_d_q(j, 1)               qair1(j) = yq(j, 1) + y_d_q(j, 1)
452               zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &               zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
453                    1)))*(ypaprs(j, 1)-ypplay(j, 1))                    1))) * (ypaprs(j, 1)-ypplay(j, 1))
454               tairsol(j) = yts(j) + y_d_ts(j)               tairsol(j) = yts(j) + y_d_ts(j)
455               rugo1(j) = yrugos(j)               rugo1(j) = yrugos(j)
456               IF (nsrf == is_oce) THEN               IF (nsrf == is_oce) THEN
457                  rugo1(j) = rugos(i, nsrf)                  rugo1(j) = frugs(i, nsrf)
458               END IF               END IF
459               psfce(j) = ypaprs(j, 1)               psfce(j) = ypaprs(j, 1)
460               patm(j) = ypplay(j, 1)               patm(j) = ypplay(j, 1)
# Line 571  contains Line 462  contains
462               qairsol(j) = yqsurf(j)               qairsol(j) = yqsurf(j)
463            END DO            END DO
464    
465            CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &            CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), qair1, &
466                 zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &                 zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, &
467                 yt10m, yq10m, yu10m, yustar)                 yq10m, wind10m(:knon), ustar(:knon))
468    
469            DO j = 1, knon            DO j = 1, knon
470               i = ni(j)               i = ni(j)
471               t2m(i, nsrf) = yt2m(j)               t2m(i, nsrf) = yt2m(j)
472               q2m(i, nsrf) = yq2m(j)               q2m(i, nsrf) = yq2m(j)
473    
474               ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman               u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
475               u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)                    / sqrt(u1(j)**2 + v1(j)**2)
476               v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)               v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
477                      / sqrt(u1(j)**2 + v1(j)**2)
478            END DO            END DO
479    
480            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, yustar, y_flux_t(:knon), &            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
481                 y_flux_q(:knon), yu, yv, yt, yq, ypblh(:knon), ycapcl, &                 y_flux_q(:knon), yu, yv, yt, yq, ypblh(:knon), ycapcl, &
482                 yoliqcl, ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)                 yoliqcl, ycteicl, ypblt, ytherm, ylcl)
483    
484            DO j = 1, knon            DO j = 1, knon
485               i = ni(j)               i = ni(j)
# Line 599  contains Line 490  contains
490               cteicl(i, nsrf) = ycteicl(j)               cteicl(i, nsrf) = ycteicl(j)
491               pblt(i, nsrf) = ypblt(j)               pblt(i, nsrf) = ypblt(j)
492               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)  
493            END DO            END DO
494    
495            DO j = 1, knon            DO j = 1, knon
# Line 610  contains Line 498  contains
498                  q2(i, k, nsrf) = yq2(j, k)                  q2(i, k, nsrf) = yq2(j, k)
499               END DO               END DO
500            END DO            END DO
501           else
502              fsnow(:, nsrf) = 0.
503         end IF if_knon         end IF if_knon
504      END DO loop_surface      END DO loop_surface
505    
506      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
507      rugos(:, is_oce) = rugmer      frugs(:, is_oce) = rugmer
508      pctsrf(:, is_oce) = pctsrf_new_oce      pctsrf(:, is_oce) = pctsrf_new_oce
509      pctsrf(:, is_sic) = pctsrf_new_sic      pctsrf(:, is_sic) = pctsrf_new_sic
510    
511      firstcal = .false.      firstcal = .false.
512    
513    END SUBROUTINE clmain    END SUBROUTINE pbl_surface
514    
515  end module clmain_m  end module pbl_surface_m

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

  ViewVC Help
Powered by ViewVC 1.1.21