/[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 186 by guez, Mon Mar 21 15:36:26 2016 UTC trunk/phylmd/Interface_surf/pbl_surface.f revision 293 by guez, Wed Jul 25 16:16:53 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, itap, pctsrf, pctsrf_new, t, q, u, v, jour, rmu0, &    SUBROUTINE pbl_surface(dtime, pctsrf, t, q, u, v, julien, mu0, ftsol, &
8         ts, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, &         cdmmax, cdhmax, ftsoil, qsol, paprs, pplay, fsnow, qsurf, evap, falbe, &
9         paprs, pplay, snow, qsurf, evap, falbe, fluxlat, rain_fall, snow_f, &         fluxlat, rain_fall, snow_f, fsolsw, fsollw, frugs, agesno, rugoro, d_t, &
10         solsw, sollw, fder, rlat, rugos, debut, agesno, rugoro, d_t, d_q, d_u, &         d_q, d_u, d_v, d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, &
11         d_v, d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, &         q2, dflux_t, dflux_q, coefh, t2m, q2m, u10m_srf, v10m_srf, pblh, capcl, &
12         dflux_t, dflux_q, ycoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh, capcl, &         oliqcl, cteicl, pblt, therm, plcl, fqcalving, ffonte, run_off_lic_0)
        oliqcl, cteicl, pblt, therm, 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 22  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  
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
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 ustarhb_m, only: ustarhb      use time_phylmdz, only: itap
     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)
     INTEGER, INTENT(IN):: itap ! numero du pas de temps  
     REAL, INTENT(inout):: pctsrf(klon, nbsrf)  
40    
41      ! la nouvelle repartition des surfaces sortie de l'interface      REAL, INTENT(inout):: pctsrf(klon, nbsrf)
42      REAL, INTENT(out):: pctsrf_new(klon, nbsrf)      ! 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 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
71    
72      REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)      REAL, INTENT(IN):: fsolsw(klon, nbsrf), fsollw(klon, nbsrf)
73      REAL, intent(in):: fder(klon)      REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)
     REAL, INTENT(IN):: rlat(klon) ! latitude en degr\'es  
   
     REAL rugos(klon, nbsrf)  
     ! rugos----input-R- longeur de rugosite (en m)  
   
     LOGICAL, INTENT(IN):: debut  
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)
86        ! flux de chaleur sensible (Cp T) (W / m2) (orientation positive vers
87        ! le bas) à la surface
88    
89        REAL, intent(out):: flux_q(klon, nbsrf)
90        ! flux de vapeur d'eau (kg / m2 / s) à la surface
91    
92      REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)      REAL, intent(out):: flux_u(klon, nbsrf), flux_v(klon, nbsrf)
93      ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)      ! tension du vent (flux turbulent de vent) à la surface, en Pa
     !                    (orientation positive vers le bas)  
     ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)  
   
     REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)  
     ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal  
     ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal  
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 : pbl, hbtm (Comme les autres  
110      ! diagnostics on cumule dans physiq ce qui permet de sortir les      REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)
111      ! grandeurs par sous-surface)      ! composantes du vent \`a 10m sans spirale d'Ekman
112      REAL pblh(klon, nbsrf)  
113      ! pblh------- HCL      ! 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
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:
133    
134        LOGICAL:: firstcal = .true.
135    
136        ! la nouvelle repartition des surfaces sortie de l'interface
137        REAL, save:: pctsrf_new_oce(klon)
138        REAL, save:: pctsrf_new_sic(klon)
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) ! longueur 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, klev), y_flux_q(klon, klev)      REAL y_flux_t(klon), y_flux_q(klon)
157      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)      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 extrapola.      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 220  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 240  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 256  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.
216      yu = 0.      yu = 0.
217      yv = 0.      yv = 0.
     yt = 0.  
218      yq = 0.      yq = 0.
     pctsrf_new = 0.  
     y_flux_u = 0.  
     y_flux_v = 0.  
     y_dflux_t = 0.  
     y_dflux_q = 0.  
     ytsoil = 999999.  
219      yrugoro = 0.      yrugoro = 0.
220      d_ts = 0.      d_ts = 0.
     yfluxlat = 0.  
221      flux_t = 0.      flux_t = 0.
222      flux_q = 0.      flux_q = 0.
223      flux_u = 0.      flux_u = 0.
224      flux_v = 0.      flux_v = 0.
225        fluxlat = 0.
226      d_t = 0.      d_t = 0.
227      d_q = 0.      d_q = 0.
228      d_u = 0.      d_u = 0.
229      d_v = 0.      d_v = 0.
230      ycoefh = 0.      coefh = 0.
231        fqcalving = 0.
232    
233      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
234      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
235      ! (\`a affiner)      ! (\`a affiner)
236    
237      pctsrf_pot = pctsrf      pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
238        pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
239      pctsrf_pot(:, is_oce) = 1. - zmasq      pctsrf_pot(:, is_oce) = 1. - zmasq
240      pctsrf_pot(:, is_sic) = 1. - zmasq      pctsrf_pot(:, is_sic) = 1. - zmasq
241    
242        ! Tester si c'est le moment de lire le fichier:
243        if (mod(itap - 1, lmt_pas) == 0) then
244           CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
245        endif
246    
247      ! Boucler sur toutes les sous-fractions du sol:      ! Boucler sur toutes les sous-fractions du sol:
248    
249      loop_surface: DO nsrf = 1, nbsrf      loop_surface: DO nsrf = 1, nbsrf
# Line 322  contains Line 263  contains
263            DO j = 1, knon            DO j = 1, knon
264               i = ni(j)               i = ni(j)
265               ypct(j) = pctsrf(i, nsrf)               ypct(j) = pctsrf(i, nsrf)
266               yts(j) = ts(i, nsrf)               yts(j) = ftsol(i, nsrf)
267               ysnow(j) = snow(i, nsrf)               snow(j) = fsnow(i, nsrf)
268               yqsurf(j) = qsurf(i, nsrf)               yqsurf(j) = qsurf(i, nsrf)
269               yalb(j) = falbe(i, nsrf)               yalb(j) = falbe(i, nsrf)
270               yrain_f(j) = rain_fall(i)               yrain_f(j) = rain_fall(i)
271               ysnow_f(j) = snow_f(i)               ysnow_f(j) = snow_f(i)
272               yagesno(j) = agesno(i, nsrf)               yagesno(j) = agesno(i, nsrf)
273               yfder(j) = fder(i)               yrugos(j) = frugs(i, nsrf)
              yrugos(j) = rugos(i, nsrf)  
274               yrugoro(j) = rugoro(i)               yrugoro(j) = rugoro(i)
275               yu1(j) = u1lay(i)               yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)
276               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)  
277               y_run_off_lic_0(j) = run_off_lic_0(i)               y_run_off_lic_0(j) = run_off_lic_0(i)
278            END DO            END DO
279    
280            ! For continent, copy soil water content            ! For continent, copy soil water content
281            IF (nsrf == is_ter) THEN            IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
              yqsol(:knon) = qsol(ni(:knon))  
           ELSE  
              yqsol = 0.  
           END IF  
282    
283            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  
284    
285            DO k = 1, klev            DO k = 1, klev
286               DO j = 1, knon               DO j = 1, knon
# Line 366  contains Line 295  contains
295               END DO               END DO
296            END DO            END DO
297    
298            ! calculer Cdrag et les coefficients d'echange            ! Calculer les géopotentiels de chaque couche:
299            CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &  
300                 yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))            zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
301                   + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
302    
303              DO k = 2, klev
304                 zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
305                      * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
306                      * (ypplay(:knon, k - 1) - ypplay(:knon, k))
307              ENDDO
308    
309              CALL cdrag(nsrf, sqrt(yu(:knon, 1)**2 + yv(:knon, 1)**2), &
310                   yt(:knon, 1), yq(:knon, 1), zgeop(:knon, 1), ypaprs(:knon, 1), &
311                   yts(:knon), yqsurf(:knon), yrugos(:knon), ycdragm(:knon), &
312                   ycdragh(:knon))
313    
314            IF (iflag_pbl == 1) THEN            IF (iflag_pbl == 1) THEN
315               CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)               ycdragm(:knon) = max(ycdragm(:knon), 0.)
316               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))               ycdragh(:knon) = max(ycdragh(:knon), 0.)
317               coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))            end IF
           END IF  
318    
319            ! on met un seuil pour coefm et coefh            ! on met un seuil pour ycdragm et ycdragh
320            IF (nsrf == is_oce) THEN            IF (nsrf == is_oce) THEN
321               coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)               ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
322               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, :))  
323            END IF            END IF
324    
325            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  
326               DO k = 1, klev + 1               DO k = 1, klev + 1
327                  DO j = 1, knon                  DO j = 1, knon
328                     i = ni(j)                     i = ni(j)
329                     yq2(j, k) = q2(i, k, nsrf)                     yq2(j, k) = q2(i, k, nsrf)
330                  END DO                  END DO
331               END DO               END DO
332              end IF
333    
334               CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)            call coef_diff_turb(dtime, nsrf, ni(:knon), ypaprs(:knon, :), &
335               IF (prt_level > 9) PRINT *, 'USTAR = ', yustar                 ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
336                   yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
337               ! iflag_pbl peut \^etre utilis\'e comme longueur de m\'elange                 ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
338    
339               IF (iflag_pbl >= 11) THEN            CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
340                  CALL vdif_kcay(knon, dtime, rg, ypaprs, yzlev, yzlay, yu, yv, &                 ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
341                       yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, yustar, &                 ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
342                       iflag_pbl)                 y_flux_u(:knon))
343               ELSE            CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
344                  CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &                 ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
345                       coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)                 ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
346               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)  
           CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &  
                ypplay, ydelp, y_d_v, y_flux_v)  
347    
348            ! calculer la diffusion de "q" et de "h"            ! calculer la diffusion de "q" et de "h"
349            CALL clqh(dtime, itap, jour, debut, rlat, knon, nsrf, ni(:knon), &            CALL clqh(dtime, julien, firstcal, nsrf, ni(:knon), &
350                 pctsrf, ytsoil, yqsol, rmu0, yrugos, yrugoro, yu1, &                 ytsoil(:knon, :), yqsol(:knon), mu0, yrugos(:knon), &
351                 yv1, coefh(:knon, :), yt, yq, yts, ypaprs, ypplay, ydelp, &                 yrugoro(:knon), yu(:knon, 1), yv(:knon, 1), ycoefh(:knon, :), &
352                 yrads, yalb(:knon), ysnow, yqsurf, yrain_f, ysnow_f, yfder, &                 ycdragh(:knon), yt(:knon, :), yq(:knon, :), yts(:knon), &
353                 yfluxlat, pctsrf_new, yagesno(:knon), y_d_t, y_d_q, &                 ypaprs(:knon, :), ypplay(:knon, :), ydelp(:knon, :), &
354                 y_d_ts(:knon), yz0_new, y_flux_t, y_flux_q, y_dflux_t, &                 yrads(:knon), yalb(:knon), snow(:knon), yqsurf(:knon), yrain_f, &
355                 y_dflux_q, y_fqcalving, y_ffonte, y_run_off_lic_0)                 ysnow_f, yfluxlat(:knon), pctsrf_new_sic, yagesno(:knon), &
356                   y_d_t(:knon, :), y_d_q(:knon, :), y_d_ts(:knon), &
357                   yz0_new(:knon), y_flux_t(:knon), y_flux_q(:knon), &
358                   y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving(:knon), &
359                   y_ffonte, y_run_off_lic_0)
360    
361            ! calculer la longueur de rugosite sur ocean            ! calculer la longueur de rugosite sur ocean
362    
363            yrugm = 0.            yrugm = 0.
364    
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
           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  
373    
374            DO k = 1, klev            DO k = 1, klev
375               DO j = 1, knon               DO j = 1, knon
376                  i = ni(j)                  i = ni(j)
377                  coefh(j, k) = coefh(j, k)*ypct(j)                  y_d_t(j, k) = y_d_t(j, k) * ypct(j)
378                  coefm(j, k) = coefm(j, k)*ypct(j)                  y_d_q(j, k) = y_d_q(j, k) * ypct(j)
379                  y_d_t(j, k) = y_d_t(j, k)*ypct(j)                  y_d_u(j, k) = y_d_u(j, k) * ypct(j)
380                  y_d_q(j, k) = y_d_q(j, k)*ypct(j)                  y_d_v(j, k) = y_d_v(j, k) * ypct(j)
                 flux_t(i, k, nsrf) = y_flux_t(j, k)  
                 flux_q(i, k, nsrf) = y_flux_q(j, k)  
                 flux_u(i, k, nsrf) = y_flux_u(j, k)  
                 flux_v(i, k, nsrf) = y_flux_v(j, k)  
                 y_d_u(j, k) = y_d_u(j, k)*ypct(j)  
                 y_d_v(j, k) = y_d_v(j, k)*ypct(j)  
381               END DO               END DO
382            END DO            END DO
383    
384            evap(:, nsrf) = -flux_q(:, 1, nsrf)            flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
385              flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
386              flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
387              flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
388    
389              evap(:, nsrf) = -flux_q(:, nsrf)
390    
391            falbe(:, nsrf) = 0.            falbe(:, nsrf) = 0.
392            snow(:, nsrf) = 0.            fsnow(:, nsrf) = 0.
393            qsurf(:, nsrf) = 0.            qsurf(:, nsrf) = 0.
394            rugos(:, nsrf) = 0.            frugs(:, nsrf) = 0.
           fluxlat(:, nsrf) = 0.  
395            DO j = 1, knon            DO j = 1, knon
396               i = ni(j)               i = ni(j)
397               d_ts(i, nsrf) = y_d_ts(j)               d_ts(i, nsrf) = y_d_ts(j)
398               falbe(i, nsrf) = yalb(j)               falbe(i, nsrf) = yalb(j)
399               snow(i, nsrf) = ysnow(j)               fsnow(i, nsrf) = snow(j)
400               qsurf(i, nsrf) = yqsurf(j)               qsurf(i, nsrf) = yqsurf(j)
401               rugos(i, nsrf) = yz0_new(j)               frugs(i, nsrf) = yz0_new(j)
402               fluxlat(i, nsrf) = yfluxlat(j)               fluxlat(i, nsrf) = yfluxlat(j)
403               IF (nsrf == is_oce) THEN               IF (nsrf == is_oce) THEN
404                  rugmer(i) = yrugm(j)                  rugmer(i) = yrugm(j)
405                  rugos(i, nsrf) = yrugm(j)                  frugs(i, nsrf) = yrugm(j)
406               END IF               END IF
407               agesno(i, nsrf) = yagesno(j)               agesno(i, nsrf) = yagesno(j)
408               fqcalving(i, nsrf) = y_fqcalving(j)               fqcalving(i, nsrf) = y_fqcalving(j)
409               ffonte(i, nsrf) = y_ffonte(j)               ffonte(i, nsrf) = y_ffonte(j)
410               cdragh(i) = cdragh(i) + coefh(j, 1)               cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
411               cdragm(i) = cdragm(i) + coefm(j, 1)               cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
412               dflux_t(i) = dflux_t(i) + y_dflux_t(j)               dflux_t(i) = dflux_t(i) + y_dflux_t(j) * ypct(j)
413               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)  
414            END DO            END DO
415            IF (nsrf == is_ter) THEN            IF (nsrf == is_ter) THEN
416               qsol(ni(:knon)) = yqsol(:knon)               qsol(ni(:knon)) = yqsol(:knon)
# Line 522  contains Line 422  contains
422            END IF            END IF
423    
424            ftsoil(:, :, nsrf) = 0.            ftsoil(:, :, nsrf) = 0.
425            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  
426    
427            DO j = 1, knon            DO j = 1, knon
428               i = ni(j)               i = ni(j)
# Line 536  contains Line 431  contains
431                  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)
432                  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)
433                  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)  
434               END DO               END DO
435            END DO            END DO
436    
437              forall (k = 2:klev) coefh(ni(:knon), k) &
438                   = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypct(:knon)
439    
440            ! diagnostic t, q a 2m et u, v a 10m            ! diagnostic t, q a 2m et u, v a 10m
441    
442            DO j = 1, knon            DO j = 1, knon
443               i = ni(j)               i = ni(j)
444               uzon(j) = yu(j, 1) + y_d_u(j, 1)               u1(j) = yu(j, 1) + y_d_u(j, 1)
445               vmer(j) = yv(j, 1) + y_d_v(j, 1)               v1(j) = yv(j, 1) + y_d_v(j, 1)
446               tair1(j) = yt(j, 1) + y_d_t(j, 1)               tair1(j) = yt(j, 1) + y_d_t(j, 1)
447               qair1(j) = yq(j, 1) + y_d_q(j, 1)               qair1(j) = yq(j, 1) + y_d_q(j, 1)
448               zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &               zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
449                    1)))*(ypaprs(j, 1)-ypplay(j, 1))                    1))) * (ypaprs(j, 1)-ypplay(j, 1))
450               tairsol(j) = yts(j) + y_d_ts(j)               tairsol(j) = yts(j) + y_d_ts(j)
451               rugo1(j) = yrugos(j)               rugo1(j) = yrugos(j)
452               IF (nsrf == is_oce) THEN               IF (nsrf == is_oce) THEN
453                  rugo1(j) = rugos(i, nsrf)                  rugo1(j) = frugs(i, nsrf)
454               END IF               END IF
455               psfce(j) = ypaprs(j, 1)               psfce(j) = ypaprs(j, 1)
456               patm(j) = ypplay(j, 1)               patm(j) = ypplay(j, 1)
# Line 561  contains Line 458  contains
458               qairsol(j) = yqsurf(j)               qairsol(j) = yqsurf(j)
459            END DO            END DO
460    
461            CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &            CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), qair1, &
462                 zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &                 zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, &
463                 yt10m, yq10m, yu10m, yustar)                 yq10m, wind10m(:knon), ustar(:knon))
464    
465            DO j = 1, knon            DO j = 1, knon
466               i = ni(j)               i = ni(j)
467               t2m(i, nsrf) = yt2m(j)               t2m(i, nsrf) = yt2m(j)
468               q2m(i, nsrf) = yq2m(j)               q2m(i, nsrf) = yq2m(j)
469    
470               ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman               u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
471               u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)                    / sqrt(u1(j)**2 + v1(j)**2)
472               v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)               v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
473                      / sqrt(u1(j)**2 + v1(j)**2)
474            END DO            END DO
475    
476            CALL hbtm(knon, ypaprs, ypplay, yt2m, yq2m, yustar, y_flux_t, &            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
477                 y_flux_q, yu, yv, yt, yq, ypblh(:knon), ycapcl, yoliqcl, &                 y_flux_q(:knon), yu, yv, yt(:knon, :), yq, ypblh(:knon), &
478                 ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)                 ycapcl, yoliqcl, ycteicl, ypblt, ytherm, ylcl)
479    
480            DO j = 1, knon            DO j = 1, knon
481               i = ni(j)               i = ni(j)
# Line 589  contains Line 486  contains
486               cteicl(i, nsrf) = ycteicl(j)               cteicl(i, nsrf) = ycteicl(j)
487               pblt(i, nsrf) = ypblt(j)               pblt(i, nsrf) = ypblt(j)
488               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)  
489            END DO            END DO
490    
491            DO j = 1, knon            DO j = 1, knon
# Line 600  contains Line 494  contains
494                  q2(i, k, nsrf) = yq2(j, k)                  q2(i, k, nsrf) = yq2(j, k)
495               END DO               END DO
496            END DO            END DO
497           else
498              fsnow(:, nsrf) = 0.
499         end IF if_knon         end IF if_knon
500      END DO loop_surface      END DO loop_surface
501    
502      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
503        frugs(:, is_oce) = rugmer
504        pctsrf(:, is_oce) = pctsrf_new_oce
505        pctsrf(:, is_sic) = pctsrf_new_sic
506    
507      rugos(:, is_oce) = rugmer      firstcal = .false.
     pctsrf = pctsrf_new  
508    
509    END SUBROUTINE clmain    END SUBROUTINE pbl_surface
510    
511  end module clmain_m  end module pbl_surface_m

Legend:
Removed from v.186  
changed lines
  Added in v.293

  ViewVC Help
Powered by ViewVC 1.1.21