/[lmdze]/trunk/phylmd/pbl_surface.f
ViewVC logotype

Diff of /trunk/phylmd/pbl_surface.f

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21