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

Legend:
Removed from v.209  
changed lines
  Added in v.298

  ViewVC Help
Powered by ViewVC 1.1.21