/[lmdze]/trunk/phylmd/Interface_surf/pbl_surface.f90
ViewVC logotype

Diff of /trunk/phylmd/Interface_surf/pbl_surface.f90

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

trunk/Sources/phylmd/clmain.f revision 145 by guez, Tue Jun 16 15:23:29 2015 UTC trunk/phylmd/Interface_surf/pbl_surface.f90 revision 343 by guez, Mon Oct 28 08:14:26 2019 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(pctsrf, t_seri, q_seri, u_seri, v_seri, julien, mu0, &
8         co2_ppm, ts, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, &         ftsol, cdmmax, cdhmax, ftsoil, qsol, paprs, play, fsnow, fqsurf, falbe, &
9         paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat, rain_fall, &         fluxlat, rain_fall, snow_fall, frugs, agesno, rugoro, d_t, d_q, d_u, &
10         snow_f, solsw, sollw, fder, rlat, rugos, debut, agesno, rugoro, d_t, &         d_v, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, dflux_t, &
11         d_q, d_u, d_v, d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, &         dflux_q, coefh, t2m, q2m, u10m_srf, v10m_srf, pblh, capcl, oliqcl, &
12         q2, dflux_t, dflux_q, ycoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh, &         cteicl, pblt, therm, plcl, fqcalving, ffonte, run_off_lic_0, albsol, &
13         capcl, oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl, &         sollw, solsw, tsol)
        fqcalving, ffonte, run_off_lic_0, flux_o, flux_g, tslab)  
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)
17        ! Date: Aug. 18th, 1993
18      ! Objet : interface de couche limite (diffusion verticale)      ! Objet : interface de couche limite (diffusion verticale)
19    
20      ! 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 22  contains
22      ! ne tient pas compte de la diff\'erentiation des sous-fractions      ! ne tient pas compte de la diff\'erentiation des sous-fractions
23      ! de sol.      ! de sol.
24    
25      ! 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.  
   
26      use clqh_m, only: clqh      use clqh_m, only: clqh
27      use clvent_m, only: clvent      use clvent_m, only: clvent
28      use coefkz_m, only: coefkz      use coef_diff_turb_m, only: coef_diff_turb
29      use coefkzmin_m, only: coefkzmin      USE conf_gcm_m, ONLY: lmt_pas
     USE conf_gcm_m, ONLY: prt_level  
30      USE conf_phys_m, ONLY: iflag_pbl      USE conf_phys_m, ONLY: iflag_pbl
31      USE dimens_m, ONLY: iim, jjm      USE dimphy, ONLY: klev, klon
     USE dimphy, ONLY: klev, klon, zmasq  
32      USE dimsoil, ONLY: nsoilmx      USE dimsoil, ONLY: nsoilmx
33      use hbtm_m, only: hbtm      use hbtm_m, only: hbtm
34        USE histwrite_phy_m, ONLY: histwrite_phy
35      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
36        USE interfoce_lim_m, ONLY: interfoce_lim
37        use phyetat0_m, only: masque
38      use stdlevvar_m, only: stdlevvar      use stdlevvar_m, only: stdlevvar
39      USE suphec_m, ONLY: rd, rg, rkappa      USE suphec_m, ONLY: rd, rg, rsigma
40      use ustarhb_m, only: ustarhb      use time_phylmdz, only: itap
     use vdif_kcay_m, only: vdif_kcay  
     use yamada4_m, only: yamada4  
   
     REAL, INTENT(IN):: dtime ! interval du temps (secondes)  
     INTEGER, INTENT(IN):: itap ! numero du pas de temps  
     REAL, INTENT(inout):: pctsrf(klon, nbsrf)  
41    
42      ! la nouvelle repartition des surfaces sortie de l'interface      REAL, INTENT(inout):: pctsrf(:, :) ! (klon, nbsrf)
43      REAL, INTENT(out):: pctsrf_new(klon, nbsrf)      ! pourcentages de surface de chaque maille
44    
45        REAL, INTENT(IN):: t_seri(:, :) ! (klon, klev) air temperature, in K
46        REAL, INTENT(IN):: q_seri(:, :) ! (klon, klev) mass fraction of water vapor
47        REAL, INTENT(IN):: u_seri(:, :), v_seri(:, :) ! (klon, klev) wind, in m s -1
48        INTEGER, INTENT(IN):: julien ! jour de l'annee en cours
49        REAL, intent(in):: mu0(klon) ! cosinus de l'angle solaire zenithal    
50    
51        REAL, INTENT(INout):: ftsol(:, :) ! (klon, nbsrf)
52        ! skin temperature of surface fraction, in K
53    
     REAL, INTENT(IN):: t(klon, klev) ! temperature (K)  
     REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg/kg)  
     REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse  
     INTEGER, INTENT(IN):: jour ! jour de l'annee en cours  
     REAL, intent(in):: rmu0(klon) ! cosinus de l'angle solaire zenithal      
     REAL, intent(in):: co2_ppm ! taux CO2 atmosphere  
     REAL, INTENT(IN):: ts(klon, nbsrf) ! temperature du sol (en Kelvin)  
54      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  
55    
56      REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)      REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
57      ! soil temperature of surface fraction      ! soil temperature of surface fraction
58    
59      REAL, INTENT(inout):: qsol(klon)      REAL, INTENT(inout):: qsol(:) ! (klon)
60      ! column-density of water in soil, in kg m-2      ! column-density of water in soil, in kg m-2
61    
62      REAL, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)      REAL, INTENT(IN):: paprs(klon, klev + 1) ! pression a intercouche (Pa)
63      REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)      REAL, INTENT(IN):: play(klon, klev) ! pression au milieu de couche (Pa)
     REAL snow(klon, nbsrf)  
     REAL qsurf(klon, nbsrf)  
     REAL evap(klon, nbsrf)  
     REAL albe(klon, nbsrf)  
     REAL alblw(klon, nbsrf)  
64    
65      REAL fluxlat(klon, nbsrf)      REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf)
66        ! column-density of mass of snow at the surface, in kg m-2
67    
68      REAL, intent(in):: rain_fall(klon)      REAL, INTENT(inout):: fqsurf(klon, nbsrf)
69      ! liquid water mass flux (kg/m2/s), positive down      REAL, intent(inout):: falbe(klon, nbsrf)
70    
71      REAL, intent(in):: snow_f(klon)      REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)
72      ! solid water mass flux (kg/m2/s), positive down      ! flux de chaleur latente, en W m-2
73    
74      REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)      REAL, intent(in):: rain_fall(klon)
75      REAL fder(klon)      ! liquid water mass flux (kg / m2 / s), positive down
     REAL, INTENT(IN):: rlat(klon) ! latitude en degr\'es  
76    
77      REAL rugos(klon, nbsrf)      REAL, intent(in):: snow_fall(klon)
78      ! rugos----input-R- longeur de rugosite (en m)      ! solid water mass flux (kg / m2 / s), positive down
79    
80      LOGICAL, INTENT(IN):: debut      REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)
81      real agesno(klon, nbsrf)      real, intent(inout):: agesno(:, :) ! (klon, nbsrf)
82      REAL, INTENT(IN):: rugoro(klon)      REAL, INTENT(IN):: rugoro(klon)
83    
84      REAL d_t(klon, klev), d_q(klon, klev)      REAL, intent(out):: d_t(:, :), d_q(:, :) ! (klon, klev)
85      ! d_t------output-R- le changement pour "t"      ! changement pour t_seri et q_seri
     ! d_q------output-R- le changement pour "q"  
86    
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_seri" et "v_seri"
89    
90      REAL, intent(out):: d_ts(klon, nbsrf) ! le changement pour "ts"      REAL, intent(out):: flux_t(klon, nbsrf)
91        ! flux de chaleur sensible (c_p T) (W / m2) (orientation positive
92        ! vers le bas) à la surface
93    
94      REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)      REAL, intent(out):: flux_q(klon, nbsrf)
95      ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)      ! flux de vapeur d'eau (kg / m2 / s) à la surface
96      !                    (orientation positive vers le bas)  
97      ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)      REAL, intent(out):: flux_u(:, :), flux_v(:, :) ! (klon, nbsrf)
98        ! tension du vent (flux turbulent de vent) à la surface, en Pa
     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  
99    
100      REAL, INTENT(out):: cdragh(klon), cdragm(klon)      REAL, INTENT(out):: cdragh(klon), cdragm(klon)
101      real q2(klon, klev+1, nbsrf)      real q2(klon, klev + 1, nbsrf)
102    
103      REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)      ! Ocean slab:
104      ! dflux_t derive du flux sensible      REAL, INTENT(out):: dflux_t(klon) ! derive du flux sensible
105      ! dflux_q derive du flux latent      REAL, INTENT(out):: dflux_q(klon) ! derive du flux latent
106      !IM "slab" ocean  
107        REAL, intent(out):: coefh(:, 2:) ! (klon, 2:klev)
108      REAL, intent(out):: ycoefh(klon, klev)      ! Pour pouvoir extraire les coefficients d'\'echange, le champ
109      REAL, intent(out):: zu1(klon)      ! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
110      REAL zv1(klon)      ! ce champ sur les quatre sous-surfaces du mod\`ele.
111      REAL t2m(klon, nbsrf), q2m(klon, nbsrf)  
112      REAL u10m(klon, nbsrf), v10m(klon, nbsrf)      REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)
113    
114      !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds      REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)
115      ! physiq ce qui permet de sortir les grdeurs par sous surface)      ! composantes du vent \`a 10m sans spirale d'Ekman
116      REAL pblh(klon, nbsrf)  
117      ! pblh------- HCL      ! Ionela Musat. Cf. Anne Mathieu : planetary boundary layer, hbtm.
118        ! Comme les autres diagnostics on cumule dans physiq ce qui permet
119        ! de sortir les grandeurs par sous-surface.
120        REAL pblh(klon, nbsrf) ! height of planetary boundary layer
121      REAL capcl(klon, nbsrf)      REAL capcl(klon, nbsrf)
122      REAL oliqcl(klon, nbsrf)      REAL oliqcl(klon, nbsrf)
123      REAL cteicl(klon, nbsrf)      REAL cteicl(klon, nbsrf)
124      REAL pblt(klon, nbsrf)      REAL, INTENT(inout):: pblt(:, :) ! (klon, nbsrf) temp\'erature au nveau HCL
     ! pblT------- T au nveau HCL  
125      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  
126      REAL plcl(klon, nbsrf)      REAL plcl(klon, nbsrf)
     REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)  
     ! 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  
     REAL run_off_lic_0(klon)  
   
     REAL flux_o(klon), flux_g(klon)  
     !IM "slab" ocean  
     ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')  
     ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')  
   
     REAL tslab(klon)  
     ! tslab-in/output-R temperature du slab ocean (en Kelvin)  
     ! uniqmnt pour slab  
127    
128      ! Local:      REAL, intent(out):: fqcalving(klon, nbsrf)
129        ! flux d'eau "perdue" par la surface et necessaire pour limiter la
130        ! hauteur de neige, en kg / m2 / s
131    
132      REAL y_flux_o(klon), y_flux_g(klon)      real ffonte(klon, nbsrf) ! flux thermique utilise pour fondre la neige
133      real ytslab(klon)      REAL, intent(inout):: run_off_lic_0(:) ! (klon)
     REAL y_fqcalving(klon), y_ffonte(klon)  
     real y_run_off_lic_0(klon)  
134    
135      REAL rugmer(klon)      REAL, intent(out):: albsol(:) ! (klon)
136        ! albedo du sol total, visible, moyen par maille
137    
138      REAL ytsoil(klon, nsoilmx)      REAL, intent(in):: sollw(:) ! (klon)
139        ! surface net downward longwave flux, in W m-2
140    
141      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      REAL, intent(in):: solsw(:) ! (klon)
142      REAL yalb(klon)      ! surface net downward shortwave flux, in W m-2
     REAL yalblw(klon)  
     REAL yu1(klon), yv1(klon)  
     ! on rajoute en output yu1 et yv1 qui sont les vents dans  
     ! la premiere couche  
     REAL ysnow(klon), yqsurf(klon), yagesno(klon)  
143    
144      real yqsol(klon)      REAL, intent(in):: tsol(:) ! (klon)
     ! column-density of water in soil, in kg m-2  
145    
146      REAL yrain_f(klon)      ! Local:
     ! liquid water mass flux (kg/m2/s), positive down  
147    
148      REAL ysnow_f(klon)      REAL d_ts(klon, nbsrf) ! variation of ftsol
149      ! solid water mass flux (kg/m2/s), positive down      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
150        REAL fsolsw(klon, nbsrf) ! flux solaire absorb\'e pour chaque sous-surface
151    
152      REAL ysollw(klon), ysolsw(klon)      ! la nouvelle repartition des surfaces sortie de l'interface
153      REAL yfder(klon)      REAL, save:: pctsrf_new_oce(klon)
154      REAL yrugm(klon), yrads(klon), yrugoro(klon)      REAL, save:: pctsrf_new_sic(klon)
155    
156        REAL y_fqcalving(klon), y_ffonte(klon)
157        real y_run_off_lic_0(klon), y_run_off_lic(klon)
158        REAL run_off_lic(klon) ! ruissellement total
159        REAL rugmer(klon)
160        REAL ytsoil(klon, nsoilmx)
161        REAL yts(klon), ypctsrf(klon), yz0_new(klon)
162        real yrugos(klon) ! longueur de rugosit\'e, en m
163        REAL yalb(klon)
164        REAL snow(klon) ! column-density of mass of snow at the surface, in kg m-2
165        real yqsurf(klon), yagesno(klon)
166        real yqsol(klon) ! column-density of water in soil, in kg m-2
167        REAL yrain_fall(klon) ! liquid water mass flux (kg / m2 / s), positive down
168        REAL ysnow_fall(klon) ! solid water mass flux (kg / m2 / s), positive down
169        REAL yrugm(klon), radsol(klon), yrugoro(klon)
170      REAL yfluxlat(klon)      REAL yfluxlat(klon)
   
171      REAL y_d_ts(klon)      REAL y_d_ts(klon)
172      REAL y_d_t(klon, klev), y_d_q(klon, klev)      REAL y_d_t(klon, klev), y_d_q(klon, klev)
173      REAL y_d_u(klon, klev), y_d_v(klon, klev)      REAL y_d_u(klon, klev), y_d_v(klon, klev)
174      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)      REAL y_flux_t(klon), y_flux_q(klon)
175      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)      REAL y_flux_u(klon), y_flux_v(klon)
176      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
177      REAL coefh(klon, klev), coefm(klon, klev)      REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
178        real ycdragh(klon), ycdragm(klon)
179      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
180      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
181      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)      REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
182        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)  
183      REAL delp(klon, klev)      REAL delp(klon, klev)
184      INTEGER i, k, nsrf      INTEGER i, k, nsrf
   
185      INTEGER ni(klon), knon, j      INTEGER ni(klon), knon, j
186    
187      REAL pctsrf_pot(klon, nbsrf)      REAL pctsrf_pot(klon, nbsrf)
188      ! "pourcentage potentiel" pour tenir compte des \'eventuelles      ! "pourcentage potentiel" pour tenir compte des \'eventuelles
189      ! apparitions ou disparitions de la glace de mer      ! apparitions ou disparitions de la glace de mer
190    
191      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.      REAL yt2m(klon), yq2m(klon), wind10m(klon)
192        REAL ustar(klon)
     REAL yt2m(klon), yq2m(klon), yu10m(klon)  
     REAL yustar(klon)  
     ! -- LOOP  
     REAL yu10mx(klon)  
     REAL yu10my(klon)  
     REAL ywindsp(klon)  
     ! -- LOOP  
193    
194      REAL yt10m(klon), yq10m(klon)      REAL yt10m(klon), yq10m(klon)
195      REAL ypblh(klon)      REAL ypblh(klon)
# Line 240  contains Line 199  contains
199      REAL ycteicl(klon)      REAL ycteicl(klon)
200      REAL ypblt(klon)      REAL ypblt(klon)
201      REAL ytherm(klon)      REAL ytherm(klon)
202      REAL ytrmb1(klon)      REAL u1(klon), v1(klon)
203      REAL ytrmb2(klon)      REAL tair1(klon)
     REAL ytrmb3(klon)  
     REAL uzon(klon), vmer(klon)  
     REAL tair1(klon), qair1(klon), tairsol(klon)  
     REAL psfce(klon), patm(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    
209        albsol = sum(falbe * pctsrf, dim = 2)
210    
211        ! R\'epartition sous maille des flux longwave et shortwave
212        ! R\'epartition du longwave par sous-surface lin\'earis\'ee
213    
214        forall (nsrf = 1:nbsrf)
215           fsollw(:, nsrf) = sollw + 4. * RSIGMA * tsol**3 &
216                * (tsol - ftsol(:, nsrf))
217           fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
218        END forall
219    
220      ytherm = 0.      ytherm = 0.
221    
222      DO k = 1, klev ! epaisseur de couche      DO k = 1, klev ! epaisseur de couche
223         DO i = 1, klon         DO i = 1, klon
224            delp(i, k) = paprs(i, k) - paprs(i, k+1)            delp(i, k) = paprs(i, k) - paprs(i, k + 1)
225         END DO         END DO
226      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  
227    
228      ! Initialization:      ! Initialization:
229      rugmer = 0.      rugmer = 0.
# Line 276  contains Line 231  contains
231      cdragm = 0.      cdragm = 0.
232      dflux_t = 0.      dflux_t = 0.
233      dflux_q = 0.      dflux_q = 0.
     zu1 = 0.  
     zv1 = 0.  
     ypct = 0.  
     yts = 0.  
     ysnow = 0.  
     yqsurf = 0.  
     yalb = 0.  
     yalblw = 0.  
     yrain_f = 0.  
     ysnow_f = 0.  
     yfder = 0.  
     ysolsw = 0.  
     ysollw = 0.  
     yrugos = 0.  
     yu1 = 0.  
     yv1 = 0.  
     yrads = 0.  
234      ypaprs = 0.      ypaprs = 0.
235      ypplay = 0.      ypplay = 0.
236      ydelp = 0.      ydelp = 0.
     yu = 0.  
     yv = 0.  
     yt = 0.  
     yq = 0.  
     pctsrf_new = 0.  
     y_flux_u = 0.  
     y_flux_v = 0.  
     y_dflux_t = 0.  
     y_dflux_q = 0.  
     ytsoil = 999999.  
237      yrugoro = 0.      yrugoro = 0.
     yu10mx = 0.  
     yu10my = 0.  
     ywindsp = 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        fqcalving = 0.
250        run_off_lic = 0.
251    
252      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
253      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
254      ! (\`a affiner)      ! (\`a affiner).
255    
256      pctsrf_pot = pctsrf      pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
257      pctsrf_pot(:, is_oce) = 1. - zmasq      pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
258      pctsrf_pot(:, is_sic) = 1. - zmasq      pctsrf_pot(:, is_oce) = 1. - masque
259        pctsrf_pot(:, is_sic) = 1. - masque
260    
261        ! Tester si c'est le moment de lire le fichier:
262        if (mod(itap - 1, lmt_pas) == 0) then
263           CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
264        endif
265    
266      ! Boucler sur toutes les sous-fractions du sol:      ! Boucler sur toutes les sous-fractions du sol:
267    
268      loop_surface: DO nsrf = 1, nbsrf      loop_surface: DO nsrf = 1, nbsrf
269         ! Chercher les indices :         ! Define ni and knon:
270    
271         ni = 0         ni = 0
272         knon = 0         knon = 0
273    
274         DO i = 1, klon         DO i = 1, klon
275            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
276            ! "potentielles"            ! "potentielles"
# Line 346  contains Line 281  contains
281         END DO         END DO
282    
283         if_knon: IF (knon /= 0) then         if_knon: IF (knon /= 0) then
284            DO j = 1, knon            ypctsrf(:knon) = pctsrf(ni(:knon), nsrf)
285               i = ni(j)            yts(:knon) = ftsol(ni(:knon), nsrf)
286               ypct(j) = pctsrf(i, nsrf)            snow(:knon) = fsnow(ni(:knon), nsrf)
287               yts(j) = ts(i, nsrf)            yqsurf(:knon) = fqsurf(ni(:knon), nsrf)
288               ytslab(i) = tslab(i)            yalb(:knon) = falbe(ni(:knon), nsrf)
289               ysnow(j) = snow(i, nsrf)            yrain_fall(:knon) = rain_fall(ni(:knon))
290               yqsurf(j) = qsurf(i, nsrf)            ysnow_fall(:knon) = snow_fall(ni(:knon))
291               yalb(j) = albe(i, nsrf)            yagesno(:knon) = agesno(ni(:knon), nsrf)
292               yalblw(j) = alblw(i, nsrf)            yrugos(:knon) = frugs(ni(:knon), nsrf)
293               yrain_f(j) = rain_fall(i)            yrugoro(:knon) = rugoro(ni(:knon))
294               ysnow_f(j) = snow_f(i)            radsol(:knon) = fsolsw(ni(:knon), nsrf) + fsollw(ni(:knon), nsrf)
295               yagesno(j) = agesno(i, nsrf)            ypaprs(:knon, klev + 1) = paprs(ni(:knon), klev + 1)
296               yfder(j) = fder(i)            y_run_off_lic_0(:knon) = run_off_lic_0(ni(:knon))
              ysolsw(j) = solsw(i, nsrf)  
              ysollw(j) = sollw(i, nsrf)  
              yrugos(j) = rugos(i, nsrf)  
              yrugoro(j) = rugoro(i)  
              yu1(j) = u1lay(i)  
              yv1(j) = v1lay(i)  
              yrads(j) = ysolsw(j) + ysollw(j)  
              ypaprs(j, klev+1) = paprs(i, klev+1)  
              y_run_off_lic_0(j) = run_off_lic_0(i)  
              yu10mx(j) = u10m(i, nsrf)  
              yu10my(j) = v10m(i, nsrf)  
              ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))  
           END DO  
297    
298            ! For continent, copy soil water content            ! For continent, copy soil water content
299            IF (nsrf == is_ter) THEN            IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
              yqsol(:knon) = qsol(ni(:knon))  
           ELSE  
              yqsol = 0.  
           END IF  
300    
301            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  
302    
303            DO k = 1, klev            DO k = 1, klev
304               DO j = 1, knon               DO j = 1, knon
305                  i = ni(j)                  i = ni(j)
306                  ypaprs(j, k) = paprs(i, k)                  ypaprs(j, k) = paprs(i, k)
307                  ypplay(j, k) = pplay(i, k)                  ypplay(j, k) = play(i, k)
308                  ydelp(j, k) = delp(i, k)                  ydelp(j, k) = delp(i, k)
309                  yu(j, k) = u(i, k)                  yu(j, k) = u_seri(i, k)
310                  yv(j, k) = v(i, k)                  yv(j, k) = v_seri(i, k)
311                  yt(j, k) = t(i, k)                  yt(j, k) = t_seri(i, k)
312                  yq(j, k) = q(i, k)                  yq(j, k) = q_seri(i, k)
313               END DO               END DO
314            END DO            END DO
315    
316            ! calculer Cdrag et les coefficients d'echange            ! Calculer les géopotentiels de chaque couche:
           CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &  
                yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))  
           IF (iflag_pbl == 1) THEN  
              CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)  
              coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))  
              coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))  
           END IF  
   
           ! on met un seuil pour coefm et coefh  
           IF (nsrf == is_oce) THEN  
              coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)  
              coefh(:knon, 1) = min(coefh(:knon, 1), 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, :))  
           END IF  
317    
318            IF (iflag_pbl >= 3) THEN            zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
319               ! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et                 + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
              ! 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  
              DO k = 1, klev + 1  
                 DO j = 1, knon  
                    i = ni(j)  
                    yq2(j, k) = q2(i, k, nsrf)  
                 END DO  
              END DO  
320    
321               CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)            DO k = 2, klev
322               IF (prt_level > 9) PRINT *, 'USTAR = ', yustar               zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
323                      * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
324                      * (ypplay(:knon, k - 1) - ypplay(:knon, k))
325              ENDDO
326    
327               ! iflag_pbl peut \^etre utilis\'e comme longueur de m\'elange            CALL cdrag(nsrf, sqrt(yu(:knon, 1)**2 + yv(:knon, 1)**2), &
328                   yt(:knon, 1), yq(:knon, 1), zgeop(:knon, 1), ypaprs(:knon, 1), &
329                   yts(:knon), yqsurf(:knon), yrugos(:knon), ycdragm(:knon), &
330                   ycdragh(:knon))
331    
332               IF (iflag_pbl >= 11) THEN            IF (iflag_pbl == 1) THEN
333                  CALL vdif_kcay(knon, dtime, rg, ypaprs, yzlev, yzlay, yu, yv, &               ycdragm(:knon) = max(ycdragm(:knon), 0.)
334                       yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, yustar, &               ycdragh(:knon) = max(ycdragh(:knon), 0.)
335                       iflag_pbl)            end IF
              ELSE  
                 CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &  
                      coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)  
              END IF  
336    
337               coefm(:knon, 2:) = ykmm(:knon, 2:klev)            IF (nsrf == is_oce) THEN
338               coefh(:knon, 2:) = ykmn(:knon, 2:klev)               ! On met un seuil pour ycdragm et ycdragh :
339                 ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
340                 ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
341            END IF            END IF
342    
343            ! calculer la diffusion des vitesses "u" et "v"            IF (iflag_pbl >= 6) yq2(:knon, :) = q2(ni(:knon), :, nsrf)
344            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &            call coef_diff_turb(nsrf, ni(:knon), ypaprs(:knon, :), &
345                 ypplay, ydelp, y_d_u, y_flux_u)                 ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
346            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &                 yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
347                 ypplay, ydelp, y_d_v, y_flux_v)                 ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
348    
349            ! calculer la diffusion de "q" et de "h"            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
350            CALL clqh(dtime, itap, jour, debut, rlat, knon, nsrf, ni(:knon), &                 ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
351                 pctsrf, ytsoil, yqsol, rmu0, co2_ppm, yrugos, yrugoro, yu1, &                 ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
352                 yv1, coefh(:knon, :), yt, yq, yts, ypaprs, ypplay, ydelp, &                 y_flux_u(:knon))
353                 yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, yfder, &            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
354                 ysolsw, yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, &                 ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
355                 y_d_ts(:knon), yz0_new, y_flux_t, y_flux_q, y_dflux_t, &                 ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
356                 y_dflux_q, y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, &                 y_flux_v(:knon))
357                 y_flux_g)  
358              CALL clqh(julien, nsrf, ni(:knon), ytsoil(:knon, :), yqsol(:knon), &
359                   mu0(ni(:knon)), yrugos(:knon), yrugoro(:knon), yu(:knon, 1), &
360                   yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), yt(:knon, :), &
361                   yq(:knon, :), yts(:knon), ypaprs(:knon, :), ypplay(:knon, :), &
362                   ydelp(:knon, :), radsol(:knon), yalb(:knon), snow(:knon), &
363                   yqsurf(:knon), yrain_fall(:knon), ysnow_fall(:knon), &
364                   yfluxlat(:knon), pctsrf_new_sic(ni(:knon)), yagesno(:knon), &
365                   y_d_t(:knon, :), y_d_q(:knon, :), y_d_ts(:knon), &
366                   yz0_new(:knon), y_flux_t(:knon), y_flux_q(:knon), &
367                   y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving(:knon), &
368                   y_ffonte(:knon), y_run_off_lic_0(:knon), y_run_off_lic(:knon))
369    
370            ! calculer la longueur de rugosite sur ocean            ! calculer la longueur de rugosite sur ocean
371    
372            yrugm = 0.            yrugm = 0.
373    
374            IF (nsrf == is_oce) THEN            IF (nsrf == is_oce) THEN
375               DO j = 1, knon               DO j = 1, knon
376                  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) &
377                       0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))                       / rg + 0.11 * 14E-6 &
378                         / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
379                  yrugm(j) = max(1.5E-05, yrugm(j))                  yrugm(j) = max(1.5E-05, yrugm(j))
380               END DO               END DO
381            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  
382    
383            DO k = 1, klev            DO k = 1, klev
384               DO j = 1, knon               DO j = 1, knon
385                  i = ni(j)                  y_d_t(j, k) = y_d_t(j, k) * ypctsrf(j)
386                  coefh(j, k) = coefh(j, k)*ypct(j)                  y_d_q(j, k) = y_d_q(j, k) * ypctsrf(j)
387                  coefm(j, k) = coefm(j, k)*ypct(j)                  y_d_u(j, k) = y_d_u(j, k) * ypctsrf(j)
388                  y_d_t(j, k) = y_d_t(j, k)*ypct(j)                  y_d_v(j, k) = y_d_v(j, k) * ypctsrf(j)
                 y_d_q(j, k) = y_d_q(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)  
389               END DO               END DO
390            END DO            END DO
391    
392            evap(:, nsrf) = -flux_q(:, 1, nsrf)            flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
393              flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
394            albe(:, nsrf) = 0.            flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
395            alblw(:, nsrf) = 0.            flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
396            snow(:, nsrf) = 0.  
397            qsurf(:, nsrf) = 0.            falbe(:, nsrf) = 0.
398            rugos(:, nsrf) = 0.            fsnow(:, nsrf) = 0.
399            fluxlat(:, nsrf) = 0.            fqsurf(:, nsrf) = 0.
400              frugs(:, nsrf) = 0.
401            DO j = 1, knon            DO j = 1, knon
402               i = ni(j)               i = ni(j)
403               d_ts(i, nsrf) = y_d_ts(j)               d_ts(i, nsrf) = y_d_ts(j)
404               albe(i, nsrf) = yalb(j)               falbe(i, nsrf) = yalb(j)
405               alblw(i, nsrf) = yalblw(j)               fsnow(i, nsrf) = snow(j)
406               snow(i, nsrf) = ysnow(j)               fqsurf(i, nsrf) = yqsurf(j)
407               qsurf(i, nsrf) = yqsurf(j)               frugs(i, nsrf) = yz0_new(j)
              rugos(i, nsrf) = yz0_new(j)  
408               fluxlat(i, nsrf) = yfluxlat(j)               fluxlat(i, nsrf) = yfluxlat(j)
409               IF (nsrf == is_oce) THEN               IF (nsrf == is_oce) THEN
410                  rugmer(i) = yrugm(j)                  rugmer(i) = yrugm(j)
411                  rugos(i, nsrf) = yrugm(j)                  frugs(i, nsrf) = yrugm(j)
412               END IF               END IF
413               agesno(i, nsrf) = yagesno(j)               agesno(i, nsrf) = yagesno(j)
414               fqcalving(i, nsrf) = y_fqcalving(j)               fqcalving(i, nsrf) = y_fqcalving(j)
415               ffonte(i, nsrf) = y_ffonte(j)               ffonte(i, nsrf) = y_ffonte(j)
416               cdragh(i) = cdragh(i) + coefh(j, 1)               cdragh(i) = cdragh(i) + ycdragh(j) * ypctsrf(j)
417               cdragm(i) = cdragm(i) + coefm(j, 1)               cdragm(i) = cdragm(i) + ycdragm(j) * ypctsrf(j)
418               dflux_t(i) = dflux_t(i) + y_dflux_t(j)               dflux_t(i) = dflux_t(i) + y_dflux_t(j) * ypctsrf(j)
419               dflux_q(i) = dflux_q(i) + y_dflux_q(j)               dflux_q(i) = dflux_q(i) + y_dflux_q(j) * ypctsrf(j)
              zu1(i) = zu1(i) + yu1(j)  
              zv1(i) = zv1(i) + yv1(j)  
420            END DO            END DO
421            IF (nsrf == is_ter) THEN            IF (nsrf == is_ter) THEN
422               qsol(ni(:knon)) = yqsol(:knon)               qsol(ni(:knon)) = yqsol(:knon)
# Line 555  contains Line 424  contains
424               DO j = 1, knon               DO j = 1, knon
425                  i = ni(j)                  i = ni(j)
426                  run_off_lic_0(i) = y_run_off_lic_0(j)                  run_off_lic_0(i) = y_run_off_lic_0(j)
427                    run_off_lic(i) = y_run_off_lic(j)
428               END DO               END DO
429            END IF            END IF
430    
431            ftsoil(:, :, nsrf) = 0.            ftsoil(:, :, nsrf) = 0.
432            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  
433    
434            DO j = 1, knon            DO j = 1, knon
435               i = ni(j)               i = ni(j)
# Line 573  contains Line 438  contains
438                  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)
439                  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)
440                  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)  
441               END DO               END DO
442            END DO            END DO
443    
444            ! diagnostic t, q a 2m et u, v a 10m            forall (k = 2:klev) coefh(ni(:knon), k) &
445                   = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypctsrf(:knon)
446    
447            DO j = 1, knon            ! Diagnostic temp\'erature, q \`a 2 m et u, v \`a 10 m:
              i = ni(j)  
              uzon(j) = yu(j, 1) + y_d_u(j, 1)  
              vmer(j) = yv(j, 1) + y_d_v(j, 1)  
              tair1(j) = yt(j, 1) + y_d_t(j, 1)  
              qair1(j) = yq(j, 1) + y_d_q(j, 1)  
              zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &  
                   1)))*(ypaprs(j, 1)-ypplay(j, 1))  
              tairsol(j) = yts(j) + y_d_ts(j)  
              rugo1(j) = yrugos(j)  
              IF (nsrf == is_oce) THEN  
                 rugo1(j) = rugos(i, nsrf)  
              END IF  
              psfce(j) = ypaprs(j, 1)  
              patm(j) = ypplay(j, 1)  
448    
449               qairsol(j) = yqsurf(j)            u1(:knon) = yu(:knon, 1) + y_d_u(:knon, 1)
450            END DO            v1(:knon) = yv(:knon, 1) + y_d_v(:knon, 1)
451              tair1(:knon) = yt(:knon, 1) + y_d_t(:knon, 1)
452    
453              IF (nsrf == is_oce) THEN
454                 rugo1(:knon) = frugs(ni(:knon), is_oce)
455              else
456                 rugo1(:knon) = yrugos(:knon)
457              END IF
458    
459            CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &            CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), &
460                 zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &                 yq(:knon, 1) + y_d_q(:knon, 1), rd * tair1(:knon) &
461                 yt10m, yq10m, yu10m, yustar)                 / (0.5 * (ypaprs(:knon, 1) + ypplay(:knon, 1))) &
462                   * (ypaprs(:knon, 1) - ypplay(:knon, 1)), &
463                   yts(:knon) + y_d_ts(:knon), yqsurf(:knon), rugo1, &
464                   ypaprs(:knon, 1), ypplay(:knon, 1), yt2m, yq2m, yt10m, yq10m, &
465                   wind10m(:knon), ustar(:knon))
466    
467            DO j = 1, knon            DO j = 1, knon
468               i = ni(j)               i = ni(j)
469               t2m(i, nsrf) = yt2m(j)               t2m(i, nsrf) = yt2m(j)
470               q2m(i, nsrf) = yq2m(j)               q2m(i, nsrf) = yq2m(j)
471    
472               ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman               u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
473               u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)                    / sqrt(u1(j)**2 + v1(j)**2)
474               v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)               v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
475                      / sqrt(u1(j)**2 + v1(j)**2)
476            END DO            END DO
477    
478            CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
479                 y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &                 y_flux_q(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
480                 ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)                 yq(:knon, :), ypblh(:knon), ycapcl, yoliqcl, ycteicl, ypblt, &
481                   ytherm, ylcl)
482    
483            DO j = 1, knon            DO j = 1, knon
484               i = ni(j)               i = ni(j)
# Line 626  contains Line 489  contains
489               cteicl(i, nsrf) = ycteicl(j)               cteicl(i, nsrf) = ycteicl(j)
490               pblt(i, nsrf) = ypblt(j)               pblt(i, nsrf) = ypblt(j)
491               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)  
492            END DO            END DO
493    
494            DO j = 1, knon            IF (iflag_pbl >= 6) q2(ni(:knon), :, nsrf) = yq2(:knon, :)
495               DO k = 1, klev + 1         else
496                  i = ni(j)            fsnow(:, nsrf) = 0.
                 q2(i, k, nsrf) = yq2(j, k)  
              END DO  
           END DO  
           !IM "slab" ocean  
           IF (nsrf == is_oce) THEN  
              DO j = 1, knon  
                 ! on projette sur la grille globale  
                 i = ni(j)  
                 IF (pctsrf_new(i, is_oce)>epsfra) THEN  
                    flux_o(i) = y_flux_o(j)  
                 ELSE  
                    flux_o(i) = 0.  
                 END IF  
              END DO  
           END IF  
   
           IF (nsrf == is_sic) THEN  
              DO j = 1, knon  
                 i = ni(j)  
                 ! On pond\`ere lorsque l'on fait le bilan au sol :  
                 IF (pctsrf_new(i, is_sic)>epsfra) THEN  
                    flux_g(i) = y_flux_g(j)  
                 ELSE  
                    flux_g(i) = 0.  
                 END IF  
              END DO  
   
           END IF  
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        frugs(:, is_oce) = rugmer
502        pctsrf(:, is_oce) = pctsrf_new_oce
503        pctsrf(:, is_sic) = pctsrf_new_sic
504    
505        CALL histwrite_phy("run_off_lic", run_off_lic)
506        ftsol = ftsol + d_ts ! update surface temperature
507        CALL histwrite_phy("dtsvdfo", d_ts(:, is_oce))
508        CALL histwrite_phy("dtsvdft", d_ts(:, is_ter))
509        CALL histwrite_phy("dtsvdfg", d_ts(:, is_lic))
510        CALL histwrite_phy("dtsvdfi", d_ts(:, is_sic))
511    
512      rugos(:, is_oce) = rugmer    END SUBROUTINE pbl_surface
     pctsrf = pctsrf_new  
   
   END SUBROUTINE clmain  
513    
514  end module clmain_m  end module pbl_surface_m

Legend:
Removed from v.145  
changed lines
  Added in v.343

  ViewVC Help
Powered by ViewVC 1.1.21