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

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

  ViewVC Help
Powered by ViewVC 1.1.21