/[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

trunk/phylmd/clmain.f revision 104 by guez, Thu Sep 4 10:05:52 2014 UTC trunk/Sources/phylmd/clmain.f revision 251 by guez, Mon Jan 8 14:12:02 2018 UTC
# Line 4  module clmain_m Line 4  module clmain_m
4    
5  contains  contains
6    
7    SUBROUTINE clmain(dtime, itap, pctsrf, pctsrf_new, t, q, u, v, jour, rmu0, &    SUBROUTINE clmain(dtime, 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, evap, falbe, fluxlat, &
9         paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat, rain_fall, &         rain_fall, snow_f, fsolsw, fsollw, frugs, agesno, rugoro, d_t, d_q, &
10         snow_f, solsw, sollw, fder, rlat, rugos, debut, agesno, rugoro, d_t, &         d_u, d_v, d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, &
11         d_q, d_u, d_v, d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, &         dflux_t, dflux_q, coefh, t2m, q2m, u10m_srf, v10m_srf, pblh, capcl, &
12         q2, dflux_t, dflux_q, ycoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh, &         oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl, fqcalving, &
13         capcl, oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl, &         ffonte, run_off_lic_0)
        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), date: 1993/08/18
# Line 19  contains Line 18  contains
18    
19      ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul      ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
20      ! de la couche limite pour les traceurs se fait avec "cltrac" et      ! de la couche limite pour les traceurs se fait avec "cltrac" et
21      ! ne tient pas compte de la différentiation des sous-fractions de      ! ne tient pas compte de la diff\'erentiation des sous-fractions
22      ! 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.  
23    
24        use clcdrag_m, only: clcdrag
25      use clqh_m, only: clqh      use clqh_m, only: clqh
26      use clvent_m, only: clvent      use clvent_m, only: clvent
27      use coefkz_m, only: coefkz      use coef_diff_turb_m, only: coef_diff_turb
28      use coefkzmin_m, only: coefkzmin      USE conf_gcm_m, ONLY: lmt_pas
     USE conf_gcm_m, ONLY: prt_level  
29      USE conf_phys_m, ONLY: iflag_pbl      USE conf_phys_m, ONLY: iflag_pbl
     USE dimens_m, ONLY: iim, jjm  
30      USE dimphy, ONLY: klev, klon, zmasq      USE dimphy, ONLY: klev, klon, zmasq
31      USE dimsoil, ONLY: nsoilmx      USE dimsoil, ONLY: nsoilmx
32      use hbtm_m, only: hbtm      use hbtm_m, only: hbtm
33      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
34        USE interfoce_lim_m, ONLY: interfoce_lim
35      use stdlevvar_m, only: stdlevvar      use stdlevvar_m, only: stdlevvar
36      USE suphec_m, ONLY: rd, rg, rkappa      USE suphec_m, ONLY: rd, rg
37      use ustarhb_m, only: ustarhb      use time_phylmdz, only: itap
     use vdif_kcay_m, only: vdif_kcay  
     use yamada4_m, only: yamada4  
38    
39      REAL, INTENT(IN):: dtime ! interval du temps (secondes)      REAL, INTENT(IN):: dtime ! interval du temps (secondes)
     INTEGER, INTENT(IN):: itap ! numero du pas de temps  
     REAL, INTENT(inout):: pctsrf(klon, nbsrf)  
40    
41      ! la nouvelle repartition des surfaces sortie de l'interface      REAL, INTENT(inout):: pctsrf(klon, nbsrf)
42      REAL, INTENT(out):: pctsrf_new(klon, nbsrf)      ! tableau des pourcentages de surface de chaque maille
43    
44      REAL, INTENT(IN):: t(klon, klev) ! temperature (K)      REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
45      REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg/kg)      REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg / kg)
46      REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse      REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
47      INTEGER, INTENT(IN):: jour ! jour de l'annee en cours      INTEGER, INTENT(IN):: julien ! jour de l'annee en cours
48      REAL, intent(in):: rmu0(klon) ! cosinus de l'angle solaire zenithal          REAL, intent(in):: mu0(klon) ! cosinus de l'angle solaire zenithal    
49      REAL, intent(in):: 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)  
50      REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh      REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
     REAL, INTENT(IN):: ksta, ksta_ter  
     LOGICAL, INTENT(IN):: ok_kzmin  
     REAL ftsoil(klon, nsoilmx, nbsrf)  
51    
52      REAL, INTENT(inout):: qsol(klon)      REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
53        ! soil temperature of surface fraction
54    
55        REAL, INTENT(inout):: qsol(:) ! (klon)
56      ! column-density of water in soil, in kg m-2      ! column-density of water in soil, in kg m-2
57    
58      REAL, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)      REAL, INTENT(IN):: paprs(klon, klev + 1) ! pression a intercouche (Pa)
59      REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)      REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
60      REAL snow(klon, nbsrf)      REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse
61      REAL qsurf(klon, nbsrf)      REAL qsurf(klon, nbsrf)
62      REAL evap(klon, nbsrf)      REAL evap(klon, nbsrf)
63      REAL albe(klon, nbsrf)      REAL, intent(inout):: falbe(klon, nbsrf)
64      REAL alblw(klon, nbsrf)      REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)
   
     REAL fluxlat(klon, nbsrf)  
65    
66      REAL, intent(in):: rain_fall(klon)      REAL, intent(in):: rain_fall(klon)
67      ! liquid water mass flux (kg/m2/s), positive down      ! liquid water mass flux (kg / m2 / s), positive down
68    
69      REAL, intent(in):: snow_f(klon)      REAL, intent(in):: snow_f(klon)
70      ! solid water mass flux (kg/m2/s), positive down      ! solid water mass flux (kg / m2 / s), positive down
   
     REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)  
     REAL fder(klon)  
     REAL, INTENT(IN):: rlat(klon) ! latitude en degrés  
   
     REAL rugos(klon, nbsrf)  
     ! rugos----input-R- longeur de rugosite (en m)  
71    
72      LOGICAL, INTENT(IN):: debut      REAL, INTENT(IN):: fsolsw(klon, nbsrf), fsollw(klon, nbsrf)
73        REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)
74      real agesno(klon, nbsrf)      real agesno(klon, nbsrf)
75      REAL, INTENT(IN):: rugoro(klon)      REAL, INTENT(IN):: rugoro(klon)
76    
# Line 100  contains Line 81  contains
81      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
82      ! changement pour "u" et "v"      ! changement pour "u" et "v"
83    
84      REAL, intent(out):: d_ts(klon, nbsrf) ! le changement pour "ts"      REAL, intent(out):: d_ts(:, :) ! (klon, nbsrf) variation of ftsol
85    
86        REAL, intent(out):: flux_t(klon, nbsrf)
87        ! flux de chaleur sensible (Cp T) (W / m2) (orientation positive vers
88        ! le bas) à la surface
89    
90        REAL, intent(out):: flux_q(klon, nbsrf)
91        ! flux de vapeur d'eau (kg / m2 / s) à la surface
92    
93      REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)      REAL, intent(out):: flux_u(klon, nbsrf), flux_v(klon, nbsrf)
94      ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)      ! tension du vent (flux turbulent de vent) à la surface, en Pa
     !                    (orientation positive vers le bas)  
     ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)  
   
     REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)  
     ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal  
     ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal  
95    
96      REAL, INTENT(out):: cdragh(klon), cdragm(klon)      REAL, INTENT(out):: cdragh(klon), cdragm(klon)
97      real q2(klon, klev+1, nbsrf)      real q2(klon, klev + 1, nbsrf)
98    
99      REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)      REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)
100      ! dflux_t derive du flux sensible      ! dflux_t derive du flux sensible
101      ! dflux_q derive du flux latent      ! dflux_q derive du flux latent
102      !IM "slab" ocean      ! IM "slab" ocean
103    
104      REAL, intent(out):: ycoefh(klon, klev)      REAL, intent(out):: coefh(:, 2:) ! (klon, 2:klev)
105      REAL, intent(out):: zu1(klon)      ! Pour pouvoir extraire les coefficients d'\'echange, le champ
106      REAL zv1(klon)      ! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
107      REAL t2m(klon, nbsrf), q2m(klon, nbsrf)      ! ce champ sur les quatre sous-surfaces du mod\`ele.
108      REAL u10m(klon, nbsrf), v10m(klon, nbsrf)  
109        REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)
110      !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds  
111      ! physiq ce qui permet de sortir les grdeurs par sous surface)      REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)
112      REAL pblh(klon, nbsrf)      ! composantes du vent \`a 10m sans spirale d'Ekman
113      ! pblh------- HCL  
114        ! Ionela Musat. Cf. Anne Mathieu : planetary boundary layer, hbtm.
115        ! Comme les autres diagnostics on cumule dans physiq ce qui permet
116        ! de sortir les grandeurs par sous-surface.
117        REAL pblh(klon, nbsrf) ! height of planetary boundary layer
118      REAL capcl(klon, nbsrf)      REAL capcl(klon, nbsrf)
119      REAL oliqcl(klon, nbsrf)      REAL oliqcl(klon, nbsrf)
120      REAL cteicl(klon, nbsrf)      REAL cteicl(klon, nbsrf)
121      REAL pblt(klon, nbsrf)      REAL, INTENT(inout):: pblt(klon, nbsrf) ! T au nveau HCL
     ! pblT------- T au nveau HCL  
122      REAL therm(klon, nbsrf)      REAL therm(klon, nbsrf)
123      REAL trmb1(klon, nbsrf)      REAL trmb1(klon, nbsrf)
124      ! trmb1-------deep_cape      ! trmb1-------deep_cape
# Line 145  contains Line 130  contains
130      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
131      ! ffonte----Flux thermique utilise pour fondre la neige      ! ffonte----Flux thermique utilise pour fondre la neige
132      ! 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
133      !           hauteur de neige, en kg/m2/s      !           hauteur de neige, en kg / m2 / s
134      REAL run_off_lic_0(klon)      REAL run_off_lic_0(klon)
135    
     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  
   
136      ! Local:      ! Local:
137    
138      REAL y_flux_o(klon), y_flux_g(klon)      LOGICAL:: firstcal = .true.
139      real ytslab(klon)  
140        ! la nouvelle repartition des surfaces sortie de l'interface
141        REAL, save:: pctsrf_new_oce(klon)
142        REAL, save:: pctsrf_new_sic(klon)
143    
144      REAL y_fqcalving(klon), y_ffonte(klon)      REAL y_fqcalving(klon), y_ffonte(klon)
145      real y_run_off_lic_0(klon)      real y_run_off_lic_0(klon)
   
146      REAL rugmer(klon)      REAL rugmer(klon)
   
147      REAL ytsoil(klon, nsoilmx)      REAL ytsoil(klon, nsoilmx)
148        REAL yts(klon), ypct(klon), yz0_new(klon)
149      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      real yrugos(klon) ! longeur de rugosite (en m)
150      REAL yalb(klon)      REAL yalb(klon)
151      REAL yalblw(klon)      REAL snow(klon), yqsurf(klon), yagesno(klon)
152      REAL yu1(klon), yv1(klon)      real yqsol(klon) ! column-density of water in soil, in kg m-2
153      ! on rajoute en output yu1 et yv1 qui sont les vents dans      REAL yrain_f(klon) ! liquid water mass flux (kg / m2 / s), positive down
154      ! la premiere couche      REAL ysnow_f(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)  
155      REAL yrugm(klon), yrads(klon), yrugoro(klon)      REAL yrugm(klon), yrads(klon), yrugoro(klon)
   
156      REAL yfluxlat(klon)      REAL yfluxlat(klon)
   
157      REAL y_d_ts(klon)      REAL y_d_ts(klon)
158      REAL y_d_t(klon, klev), y_d_q(klon, klev)      REAL y_d_t(klon, klev), y_d_q(klon, klev)
159      REAL y_d_u(klon, klev), y_d_v(klon, klev)      REAL y_d_u(klon, klev), y_d_v(klon, klev)
160      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)      REAL y_flux_t(klon), y_flux_q(klon)
161      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)      REAL y_flux_u(klon), y_flux_v(klon)
162      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
163      REAL coefh(klon, klev), coefm(klon, klev)      REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
164        real ycdragh(klon), ycdragm(klon)
165      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
166      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
167      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)      REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
168        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)  
169      REAL delp(klon, klev)      REAL delp(klon, klev)
170      INTEGER i, k, nsrf      INTEGER i, k, nsrf
   
171      INTEGER ni(klon), knon, j      INTEGER ni(klon), knon, j
172    
173      REAL pctsrf_pot(klon, nbsrf)      REAL pctsrf_pot(klon, nbsrf)
174      ! "pourcentage potentiel" pour tenir compte des éventuelles      ! "pourcentage potentiel" pour tenir compte des \'eventuelles
175      ! apparitions ou disparitions de la glace de mer      ! apparitions ou disparitions de la glace de mer
176    
177      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.      REAL yt2m(klon), yq2m(klon), wind10m(klon)
178        REAL ustar(klon)
     REAL yt2m(klon), yq2m(klon), yu10m(klon)  
     REAL yustar(klon)  
     ! -- LOOP  
     REAL yu10mx(klon)  
     REAL yu10my(klon)  
     REAL ywindsp(klon)  
     ! -- LOOP  
179    
180      REAL yt10m(klon), yq10m(klon)      REAL yt10m(klon), yq10m(klon)
181      REAL ypblh(klon)      REAL ypblh(klon)
# Line 241  contains Line 188  contains
188      REAL ytrmb1(klon)      REAL ytrmb1(klon)
189      REAL ytrmb2(klon)      REAL ytrmb2(klon)
190      REAL ytrmb3(klon)      REAL ytrmb3(klon)
191      REAL uzon(klon), vmer(klon)      REAL u1(klon), v1(klon)
192      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
193      REAL psfce(klon), patm(klon)      REAL psfce(klon), patm(klon)
194    
195      REAL qairsol(klon), 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    
# Line 258  contains Line 202  contains
202    
203      DO k = 1, klev ! epaisseur de couche      DO k = 1, klev ! epaisseur de couche
204         DO i = 1, klon         DO i = 1, klon
205            delp(i, k) = paprs(i, k) - paprs(i, k+1)            delp(i, k) = paprs(i, k) - paprs(i, k + 1)
206         END DO         END DO
207      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  
208    
209      ! Initialization:      ! Initialization:
210      rugmer = 0.      rugmer = 0.
# Line 274  contains Line 212  contains
212      cdragm = 0.      cdragm = 0.
213      dflux_t = 0.      dflux_t = 0.
214      dflux_q = 0.      dflux_q = 0.
     zu1 = 0.  
     zv1 = 0.  
215      ypct = 0.      ypct = 0.
     yts = 0.  
     ysnow = 0.  
216      yqsurf = 0.      yqsurf = 0.
     yalb = 0.  
     yalblw = 0.  
217      yrain_f = 0.      yrain_f = 0.
218      ysnow_f = 0.      ysnow_f = 0.
     yfder = 0.  
     ysolsw = 0.  
     ysollw = 0.  
219      yrugos = 0.      yrugos = 0.
     yu1 = 0.  
     yv1 = 0.  
     yrads = 0.  
220      ypaprs = 0.      ypaprs = 0.
221      ypplay = 0.      ypplay = 0.
222      ydelp = 0.      ydelp = 0.
# Line 298  contains Line 224  contains
224      yv = 0.      yv = 0.
225      yt = 0.      yt = 0.
226      yq = 0.      yq = 0.
     pctsrf_new = 0.  
     y_flux_u = 0.  
     y_flux_v = 0.  
227      y_dflux_t = 0.      y_dflux_t = 0.
228      y_dflux_q = 0.      y_dflux_q = 0.
     ytsoil = 999999.  
229      yrugoro = 0.      yrugoro = 0.
     yu10mx = 0.  
     yu10my = 0.  
     ywindsp = 0.  
230      d_ts = 0.      d_ts = 0.
     yfluxlat = 0.  
231      flux_t = 0.      flux_t = 0.
232      flux_q = 0.      flux_q = 0.
233      flux_u = 0.      flux_u = 0.
234      flux_v = 0.      flux_v = 0.
235        fluxlat = 0.
236      d_t = 0.      d_t = 0.
237      d_q = 0.      d_q = 0.
238      d_u = 0.      d_u = 0.
239      d_v = 0.      d_v = 0.
240      ycoefh = 0.      coefh = 0.
241    
242      ! Initialisation des "pourcentages potentiels". On considère ici qu'on      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
243      ! peut avoir potentiellement de la glace sur tout le domaine océanique      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
244      ! (à affiner)      ! (\`a affiner)
245    
246      pctsrf_pot = pctsrf      pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
247        pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
248      pctsrf_pot(:, is_oce) = 1. - zmasq      pctsrf_pot(:, is_oce) = 1. - zmasq
249      pctsrf_pot(:, is_sic) = 1. - zmasq      pctsrf_pot(:, is_sic) = 1. - zmasq
250    
251        ! Tester si c'est le moment de lire le fichier:
252        if (mod(itap - 1, lmt_pas) == 0) then
253           CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
254        endif
255    
256      ! Boucler sur toutes les sous-fractions du sol:      ! Boucler sur toutes les sous-fractions du sol:
257    
258      loop_surface: DO nsrf = 1, nbsrf      loop_surface: DO nsrf = 1, nbsrf
# Line 335  contains Line 260  contains
260         ni = 0         ni = 0
261         knon = 0         knon = 0
262         DO i = 1, klon         DO i = 1, klon
263            ! Pour déterminer le domaine à traiter, on utilise les surfaces            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
264            ! "potentielles"            ! "potentielles"
265            IF (pctsrf_pot(i, nsrf) > epsfra) THEN            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
266               knon = knon + 1               knon = knon + 1
# Line 347  contains Line 272  contains
272            DO j = 1, knon            DO j = 1, knon
273               i = ni(j)               i = ni(j)
274               ypct(j) = pctsrf(i, nsrf)               ypct(j) = pctsrf(i, nsrf)
275               yts(j) = ts(i, nsrf)               yts(j) = ftsol(i, nsrf)
276               ytslab(i) = tslab(i)               snow(j) = fsnow(i, nsrf)
              ysnow(j) = snow(i, nsrf)  
277               yqsurf(j) = qsurf(i, nsrf)               yqsurf(j) = qsurf(i, nsrf)
278               yalb(j) = albe(i, nsrf)               yalb(j) = falbe(i, nsrf)
              yalblw(j) = alblw(i, nsrf)  
279               yrain_f(j) = rain_fall(i)               yrain_f(j) = rain_fall(i)
280               ysnow_f(j) = snow_f(i)               ysnow_f(j) = snow_f(i)
281               yagesno(j) = agesno(i, nsrf)               yagesno(j) = agesno(i, nsrf)
282               yfder(j) = fder(i)               yrugos(j) = frugs(i, nsrf)
              ysolsw(j) = solsw(i, nsrf)  
              ysollw(j) = sollw(i, nsrf)  
              yrugos(j) = rugos(i, nsrf)  
283               yrugoro(j) = rugoro(i)               yrugoro(j) = rugoro(i)
284               yu1(j) = u1lay(i)               yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)
285               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)  
286               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))  
287            END DO            END DO
288    
289            ! For continent, copy soil water content            ! For continent, copy soil water content
290            IF (nsrf == is_ter) THEN            IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
              yqsol(:knon) = qsol(ni(:knon))  
           ELSE  
              yqsol = 0.  
           END IF  
291    
292            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  
293    
294            DO k = 1, klev            DO k = 1, klev
295               DO j = 1, knon               DO j = 1, knon
# Line 398  contains Line 304  contains
304               END DO               END DO
305            END DO            END DO
306    
307            ! calculer Cdrag et les coefficients d'echange            ! Calculer les géopotentiels de chaque couche:
308            CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &  
309                 yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))            zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
310                   + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
311    
312              DO k = 2, klev
313                 zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
314                      * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
315                      * (ypplay(:knon, k - 1) - ypplay(:knon, k))
316              ENDDO
317    
318              CALL clcdrag(nsrf, yu(:knon, 1), yv(:knon, 1), yt(:knon, 1), &
319                   yq(:knon, 1), zgeop(:knon, 1), yts(:knon), yqsurf(:knon), &
320                   yrugos(:knon), ycdragm(:knon), ycdragh(:knon))
321    
322            IF (iflag_pbl == 1) THEN            IF (iflag_pbl == 1) THEN
323               CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)               ycdragm(:knon) = max(ycdragm(:knon), 0.)
324               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))               ycdragh(:knon) = max(ycdragh(:knon), 0.)
325               coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))            end IF
           END IF  
326    
327            ! on met un seuil pour coefm et coefh            ! on met un seuil pour ycdragm et ycdragh
328            IF (nsrf == is_oce) THEN            IF (nsrf == is_oce) THEN
329               coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)               ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
330               coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)               ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
331            END IF            END IF
332    
333            IF (ok_kzmin) THEN            IF (iflag_pbl >= 6) 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  
   
           IF (iflag_pbl >= 3) THEN  
              ! Mellor et Yamada adapté à Mars, Richard Fournier et  
              ! Frédéric 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  
334               DO k = 1, klev + 1               DO k = 1, klev + 1
335                  DO j = 1, knon                  DO j = 1, knon
336                     i = ni(j)                     i = ni(j)
337                     yq2(j, k) = q2(i, k, nsrf)                     yq2(j, k) = q2(i, k, nsrf)
338                  END DO                  END DO
339               END DO               END DO
340              end IF
341    
342               CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)            call coef_diff_turb(dtime, nsrf, ni(:knon), ypaprs(:knon, :), &
343               IF (prt_level > 9) PRINT *, 'USTAR = ', yustar                 ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
344                   yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
345               ! iflag_pbl peut être utilisé comme longueur de mélange                 ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
346    
347               IF (iflag_pbl >= 11) THEN            CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
348                  CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &                 ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
349                       yu, yv, yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, &                 ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
350                       yustar, iflag_pbl)                 y_flux_u(:knon))
351               ELSE            CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
352                  CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &                 ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
353                       coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)                 ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
354               END IF                 y_flux_v(:knon))
   
              coefm(:knon, 2:) = ykmm(:knon, 2:klev)  
              coefh(:knon, 2:) = ykmn(:knon, 2:klev)  
           END IF  
   
           ! calculer la diffusion des vitesses "u" et "v"  
           CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &  
                ypplay, ydelp, y_d_u, y_flux_u)  
           CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &  
                ypplay, ydelp, y_d_v, y_flux_v)  
355    
356            ! calculer la diffusion de "q" et de "h"            ! calculer la diffusion de "q" et de "h"
357            CALL clqh(dtime, itap, jour, debut, rlat, knon, nsrf, ni, pctsrf, &            CALL clqh(dtime, julien, firstcal, nsrf, ni(:knon), &
358                 ytsoil, yqsol, rmu0, co2_ppm, yrugos, yrugoro, &                 ytsoil(:knon, :), yqsol(:knon), mu0, yrugos, yrugoro, &
359                 yu1, yv1, coefh(:knon, :), yt, yq, yts, ypaprs, ypplay, ydelp, &                 yu(:knon, 1), yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), &
360                 yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, yfder, &                 yt, yq, yts(:knon), ypaprs, ypplay, ydelp, yrads(:knon), &
361                 ysolsw, yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts, &                 yalb(:knon), snow(:knon), yqsurf, yrain_f, ysnow_f, &
362                 yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q, &                 yfluxlat(:knon), pctsrf_new_sic, yagesno(:knon), y_d_t, y_d_q, &
363                 y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g)                 y_d_ts(:knon), yz0_new, y_flux_t(:knon), y_flux_q(:knon), &
364                   y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving, y_ffonte, &
365                   y_run_off_lic_0)
366    
367            ! calculer la longueur de rugosite sur ocean            ! calculer la longueur de rugosite sur ocean
368            yrugm = 0.            yrugm = 0.
369            IF (nsrf == is_oce) THEN            IF (nsrf == is_oce) THEN
370               DO j = 1, knon               DO j = 1, knon
371                  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) &
372                       0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))                       / rg + 0.11 * 14E-6 &
373                         / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
374                  yrugm(j) = max(1.5E-05, yrugm(j))                  yrugm(j) = max(1.5E-05, yrugm(j))
375               END DO               END DO
376            END IF            END IF
377            DO j = 1, knon            DO j = 1, knon
378               y_dflux_t(j) = y_dflux_t(j)*ypct(j)               y_dflux_t(j) = y_dflux_t(j) * ypct(j)
379               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)  
380            END DO            END DO
381    
382            DO k = 1, klev            DO k = 1, klev
383               DO j = 1, knon               DO j = 1, knon
384                  i = ni(j)                  i = ni(j)
385                  coefh(j, k) = coefh(j, k)*ypct(j)                  y_d_t(j, k) = y_d_t(j, k) * ypct(j)
386                  coefm(j, k) = coefm(j, k)*ypct(j)                  y_d_q(j, k) = y_d_q(j, k) * ypct(j)
387                  y_d_t(j, k) = y_d_t(j, k)*ypct(j)                  y_d_u(j, k) = y_d_u(j, k) * ypct(j)
388                  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)  
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              flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
395              flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
396    
397              evap(:, nsrf) = -flux_q(:, nsrf)
398    
399            albe(:, nsrf) = 0.            falbe(:, nsrf) = 0.
400            alblw(:, nsrf) = 0.            fsnow(:, nsrf) = 0.
           snow(:, nsrf) = 0.  
401            qsurf(:, nsrf) = 0.            qsurf(:, nsrf) = 0.
402            rugos(:, nsrf) = 0.            frugs(:, nsrf) = 0.
           fluxlat(:, nsrf) = 0.  
403            DO j = 1, knon            DO j = 1, knon
404               i = ni(j)               i = ni(j)
405               d_ts(i, nsrf) = y_d_ts(j)               d_ts(i, nsrf) = y_d_ts(j)
406               albe(i, nsrf) = yalb(j)               falbe(i, nsrf) = yalb(j)
407               alblw(i, nsrf) = yalblw(j)               fsnow(i, nsrf) = snow(j)
              snow(i, nsrf) = ysnow(j)  
408               qsurf(i, nsrf) = yqsurf(j)               qsurf(i, nsrf) = yqsurf(j)
409               rugos(i, nsrf) = yz0_new(j)               frugs(i, nsrf) = yz0_new(j)
410               fluxlat(i, nsrf) = yfluxlat(j)               fluxlat(i, nsrf) = yfluxlat(j)
411               IF (nsrf == is_oce) THEN               IF (nsrf == is_oce) THEN
412                  rugmer(i) = yrugm(j)                  rugmer(i) = yrugm(j)
413                  rugos(i, nsrf) = yrugm(j)                  frugs(i, nsrf) = yrugm(j)
414               END IF               END IF
415               agesno(i, nsrf) = yagesno(j)               agesno(i, nsrf) = yagesno(j)
416               fqcalving(i, nsrf) = y_fqcalving(j)               fqcalving(i, nsrf) = y_fqcalving(j)
417               ffonte(i, nsrf) = y_ffonte(j)               ffonte(i, nsrf) = y_ffonte(j)
418               cdragh(i) = cdragh(i) + coefh(j, 1)               cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
419               cdragm(i) = cdragm(i) + coefm(j, 1)               cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
420               dflux_t(i) = dflux_t(i) + y_dflux_t(j)               dflux_t(i) = dflux_t(i) + y_dflux_t(j)
421               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)  
422            END DO            END DO
423            IF (nsrf == is_ter) THEN            IF (nsrf == is_ter) THEN
424               qsol(ni(:knon)) = yqsol(:knon)               qsol(ni(:knon)) = yqsol(:knon)
# Line 554  contains Line 428  contains
428                  run_off_lic_0(i) = y_run_off_lic_0(j)                  run_off_lic_0(i) = y_run_off_lic_0(j)
429               END DO               END DO
430            END IF            END IF
431            !$$$ PB ajout pour soil  
432            ftsoil(:, :, nsrf) = 0.            ftsoil(:, :, nsrf) = 0.
433            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  
434    
435            DO j = 1, knon            DO j = 1, knon
436               i = ni(j)               i = ni(j)
# Line 570  contains Line 439  contains
439                  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)
440                  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)
441                  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)  
442               END DO               END DO
443            END DO            END DO
444    
445              forall (k = 2:klev) coefh(ni(:knon), k) &
446                   = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypct(:knon)
447    
448            ! diagnostic t, q a 2m et u, v a 10m            ! diagnostic t, q a 2m et u, v a 10m
449    
450            DO j = 1, knon            DO j = 1, knon
451               i = ni(j)               i = ni(j)
452               uzon(j) = yu(j, 1) + y_d_u(j, 1)               u1(j) = yu(j, 1) + y_d_u(j, 1)
453               vmer(j) = yv(j, 1) + y_d_v(j, 1)               v1(j) = yv(j, 1) + y_d_v(j, 1)
454               tair1(j) = yt(j, 1) + y_d_t(j, 1)               tair1(j) = yt(j, 1) + y_d_t(j, 1)
455               qair1(j) = yq(j, 1) + y_d_q(j, 1)               qair1(j) = yq(j, 1) + y_d_q(j, 1)
456               zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &               zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
457                    1)))*(ypaprs(j, 1)-ypplay(j, 1))                    1))) * (ypaprs(j, 1)-ypplay(j, 1))
458               tairsol(j) = yts(j) + y_d_ts(j)               tairsol(j) = yts(j) + y_d_ts(j)
459               rugo1(j) = yrugos(j)               rugo1(j) = yrugos(j)
460               IF (nsrf == is_oce) THEN               IF (nsrf == is_oce) THEN
461                  rugo1(j) = rugos(i, nsrf)                  rugo1(j) = frugs(i, nsrf)
462               END IF               END IF
463               psfce(j) = ypaprs(j, 1)               psfce(j) = ypaprs(j, 1)
464               patm(j) = ypplay(j, 1)               patm(j) = ypplay(j, 1)
# Line 595  contains Line 466  contains
466               qairsol(j) = yqsurf(j)               qairsol(j) = yqsurf(j)
467            END DO            END DO
468    
469            CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &            CALL stdlevvar(klon, knon, nsrf, u1(:knon), v1(:knon), tair1(:knon), &
470                 zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &                 qair1, zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, &
471                 yt10m, yq10m, yu10m, yustar)                 yq2m, yt10m, yq10m, wind10m(:knon), ustar(:knon))
472    
473            DO j = 1, knon            DO j = 1, knon
474               i = ni(j)               i = ni(j)
475               t2m(i, nsrf) = yt2m(j)               t2m(i, nsrf) = yt2m(j)
476               q2m(i, nsrf) = yq2m(j)               q2m(i, nsrf) = yq2m(j)
477    
478               ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman               u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
479               u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)                    / sqrt(u1(j)**2 + v1(j)**2)
480               v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)               v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
481                      / sqrt(u1(j)**2 + v1(j)**2)
482            END DO            END DO
483    
484            CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
485                 y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &                 y_flux_q(:knon), yu, yv, yt, yq, ypblh(:knon), ycapcl, &
486                 ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)                 yoliqcl, ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
487    
488            DO j = 1, knon            DO j = 1, knon
489               i = ni(j)               i = ni(j)
# Line 634  contains Line 505  contains
505                  q2(i, k, nsrf) = yq2(j, k)                  q2(i, k, nsrf) = yq2(j, k)
506               END DO               END DO
507            END DO            END DO
508            !IM "slab" ocean         else
509            IF (nsrf == is_oce) THEN            fsnow(:, nsrf) = 0.
              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  
510         end IF if_knon         end IF if_knon
511      END DO loop_surface      END DO loop_surface
512    
513      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
514        frugs(:, is_oce) = rugmer
515        pctsrf(:, is_oce) = pctsrf_new_oce
516        pctsrf(:, is_sic) = pctsrf_new_sic
517    
518      rugos(:, is_oce) = rugmer      firstcal = .false.
     pctsrf = pctsrf_new  
519    
520    END SUBROUTINE clmain    END SUBROUTINE clmain
521    
Legend:
Removed from v.104  
changed lines
  Added in v.251
  ViewVC Help
Powered by ViewVC 1.1.21