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

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

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

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

Legend:
Removed from v.134  
changed lines
  Added in v.307

  ViewVC Help
Powered by ViewVC 1.1.21