/[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/phylmd/clmain.f90 revision 76 by guez, Fri Nov 15 18:45:49 2013 UTC trunk/phylmd/Interface_surf/pbl_surface.f90 revision 328 by guez, Thu Jun 13 14:40:06 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, &    SUBROUTINE pbl_surface(pctsrf, t, q, u, v, julien, mu0, ftsol, cdmmax, &
8         jour, rmu0, co2_ppm, ok_veget, ocean, ts, &         cdhmax, ftsoil, qsol, paprs, play, fsnow, fqsurf, falbe, fluxlat, &
9         soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, &         rain_fall, snow_fall, frugs, agesno, rugoro, d_t, d_q, d_u, d_v, &
10         qsol, paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat, &         flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, dflux_t, dflux_q, &
11         rain_fall, snow_f, solsw, sollw, fder, rlon, rlat, &         coefh, t2m, q2m, u10m_srf, v10m_srf, pblh, capcl, oliqcl, cteicl, pblt, &
12         rugos, debut, agesno, rugoro, d_t, d_q, d_u, d_v, &         therm, plcl, fqcalving, ffonte, run_off_lic_0, albsol, sollw, solsw, &
13         d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, &         tsol)
14         dflux_t, dflux_q, ycoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh, &  
15         capcl, oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl, &      ! From phylmd/clmain.F, version 1.6, 2005/11/16 14:47:19
16         fqcalving, ffonte, run_off_lic_0, flux_o, flux_g, tslab, seaice)      ! Author: Z. X. Li (LMD/CNRS)
17        ! Date: Aug. 18th, 1993
     ! From phylmd/clmain.F, version 1.6 2005/11/16 14:47:19  
     ! Author: Z. X. Li (LMD/CNRS), date: 1993/08/18  
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 calendar, ONLY: ymds2ju      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
     USE dynetat0_m, ONLY: day_ini  
     USE gath_cpl, ONLY: gath2cpl  
33      use hbtm_m, only: hbtm      use hbtm_m, only: hbtm
34      USE histbeg_totreg_m, ONLY: histbeg_totreg      USE histwrite_phy_m, ONLY: histwrite_phy
     USE histdef_m, ONLY: histdef  
     USE histend_m, ONLY: histend  
     USE histsync_m, ONLY: histsync  
     use histwrite_m, only: histwrite  
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 suphec_m, ONLY: rd, rg, rkappa      USE interfoce_lim_m, ONLY: interfoce_lim
37      USE temps, ONLY: annee_ref, itau_phy      use phyetat0_m, only: masque
38      use ustarhb_m, only: ustarhb      use stdlevvar_m, only: stdlevvar
39      use vdif_kcay_m, only: vdif_kcay      USE suphec_m, ONLY: rd, rg, rsigma
40      use yamada4_m, only: yamada4      use time_phylmdz, only: itap
   
     ! Arguments:  
   
     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(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 co2_ppm ! taux CO2 atmosphere  
51      LOGICAL ok_veget      REAL, INTENT(INout):: ftsol(:, :) ! (klon, nbsrf)
52      CHARACTER(len=*), INTENT(IN):: ocean      ! skin temperature of surface fraction, in K
53      REAL ts(klon, nbsrf) ! input-R- temperature du sol (en Kelvin)  
     LOGICAL, INTENT(IN):: soil_model  
54      REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh      REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
     REAL ksta, ksta_ter  
     LOGICAL ok_kzmin  
     REAL ftsoil(klon, nsoilmx, nbsrf)  
     REAL qsol(klon)  
     REAL, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)  
     REAL, INTENT(IN):: pplay(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)  
   
     REAL fluxlat(klon, nbsrf)  
   
     REAL, intent(in):: rain_fall(klon), snow_f(klon)  
     REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)  
     REAL fder(klon)  
     REAL, INTENT(IN):: rlon(klon)  
     REAL, INTENT(IN):: rlat(klon) ! latitude en degrés  
55    
56      REAL rugos(klon, nbsrf)      REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
57      ! rugos----input-R- longeur de rugosite (en m)      ! soil temperature of surface fraction
58    
59        REAL, INTENT(inout):: qsol(:) ! (klon)
60        ! column-density of water in soil, in kg m-2
61    
62        REAL, INTENT(IN):: paprs(klon, klev + 1) ! pression a intercouche (Pa)
63        REAL, INTENT(IN):: play(klon, klev) ! pression au milieu de couche (Pa)
64        REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse
65        REAL, INTENT(inout):: fqsurf(klon, nbsrf)
66        REAL, intent(inout):: falbe(klon, nbsrf)
67    
68      LOGICAL, INTENT(IN):: debut      REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)
69        ! flux de chaleur latente, en W m-2
70    
71        REAL, intent(in):: rain_fall(klon)
72        ! liquid water mass flux (kg / m2 / s), positive down
73    
74        REAL, intent(in):: snow_fall(klon)
75        ! solid water mass flux (kg / m2 / s), positive down
76    
77        REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)
78      real agesno(klon, nbsrf)      real agesno(klon, nbsrf)
79      REAL, INTENT(IN):: rugoro(klon)      REAL, INTENT(IN):: rugoro(klon)
80    
81      REAL d_t(klon, klev), d_q(klon, klev)      REAL, intent(out):: d_t(:, :), d_q(:, :) ! (klon, klev)
82      ! d_t------output-R- le changement pour "t"      ! changement pour t et q
     ! d_q------output-R- le changement pour "q"  
83    
84      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
85      ! changement pour "u" et "v"      ! changement pour "u" et "v"
86    
87      REAL d_ts(klon, nbsrf)      REAL, intent(out):: flux_t(klon, nbsrf)
88      ! d_ts-----output-R- le changement pour "ts"      ! flux de chaleur sensible (c_p T) (W / m2) (orientation positive
89        ! vers le bas) à la surface
90    
91        REAL, intent(out):: flux_q(klon, nbsrf)
92        ! flux de vapeur d'eau (kg / m2 / s) à la surface
93    
94      REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)      REAL, intent(out):: flux_u(:, :), flux_v(:, :) ! (klon, nbsrf)
95      ! 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  
96    
97      REAL, INTENT(out):: cdragh(klon), cdragm(klon)      REAL, INTENT(out):: cdragh(klon), cdragm(klon)
98      real q2(klon, klev+1, nbsrf)      real q2(klon, klev + 1, nbsrf)
99    
100      REAL dflux_t(klon), dflux_q(klon)      ! Ocean slab:
101      ! dflux_t derive du flux sensible      REAL, INTENT(out):: dflux_t(klon) ! derive du flux sensible
102      ! dflux_q derive du flux latent      REAL, INTENT(out):: dflux_q(klon) ! derive du flux latent
103      !IM "slab" ocean  
104        REAL, intent(out):: coefh(:, 2:) ! (klon, 2:klev)
105      REAL, intent(out):: ycoefh(klon, klev)      ! Pour pouvoir extraire les coefficients d'\'echange, le champ
106      REAL, intent(out):: zu1(klon)      ! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
107      REAL zv1(klon)      ! ce champ sur les quatre sous-surfaces du mod\`ele.
108      REAL t2m(klon, nbsrf), q2m(klon, nbsrf)  
109      REAL u10m(klon, nbsrf), v10m(klon, nbsrf)      REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)
110    
111      !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds      REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)
112      ! physiq ce qui permet de sortir les grdeurs par sous surface)      ! composantes du vent \`a 10m sans spirale d'Ekman
113      REAL pblh(klon, nbsrf)  
114      ! pblh------- HCL      ! 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)
     REAL trmb1(klon, nbsrf)  
     ! trmb1-------deep_cape  
     REAL trmb2(klon, nbsrf)  
     ! trmb2--------inhibition  
     REAL trmb3(klon, nbsrf)  
     ! trmb3-------Point Omega  
123      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  
124    
125      REAL seaice(klon)      REAL, intent(out):: fqcalving(klon, nbsrf)
126      ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')      ! flux d'eau "perdue" par la surface et necessaire pour limiter la
127        ! hauteur de neige, en kg / m2 / s
128    
129        real ffonte(klon, nbsrf) ! flux thermique utilise pour fondre la neige
130        REAL, intent(inout):: run_off_lic_0(:) ! (klon)
131    
132        REAL, intent(out):: albsol(:) ! (klon)
133        ! albedo du sol total, visible, moyen par maille
134    
135        REAL, intent(in):: sollw(:) ! (klon)
136        ! surface net downward longwave flux, in W m-2
137    
138        REAL, intent(in):: solsw(:) ! (klon)
139        ! surface net downward shortwave flux, in W m-2
140    
141        REAL, intent(in):: tsol(:) ! (klon)
142    
143      ! Local:      ! Local:
144    
145      REAL y_flux_o(klon), y_flux_g(klon)      REAL d_ts(klon, nbsrf) ! variation of ftsol
146      real ytslab(klon)      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
147      real y_seaice(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)  
148    
149      REAL rugmer(klon)      ! la nouvelle repartition des surfaces sortie de l'interface
150        REAL, save:: pctsrf_new_oce(klon)
151        REAL, save:: pctsrf_new_sic(klon)
152    
153        REAL y_fqcalving(klon), y_ffonte(klon)
154        real y_run_off_lic_0(klon), y_run_off_lic(klon)
155        REAL run_off_lic(klon) ! ruissellement total
156        REAL rugmer(klon)
157      REAL ytsoil(klon, nsoilmx)      REAL ytsoil(klon, nsoilmx)
158        REAL yts(klon), ypctsrf(klon), yz0_new(klon)
159      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      real yrugos(klon) ! longueur de rugosite (en m)
160      REAL yalb(klon)      REAL yalb(klon)
161      REAL yalblw(klon)      REAL snow(klon), yqsurf(klon), yagesno(klon)
162      REAL yu1(klon), yv1(klon)      real yqsol(klon) ! column-density of water in soil, in kg m-2
163      ! 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
164      ! la premiere couche      REAL ysnow_fall(klon) ! solid water mass flux (kg / m2 / s), positive down
165      REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)      REAL yrugm(klon), radsol(klon), yrugoro(klon)
     REAL yrain_f(klon), ysnow_f(klon)  
     REAL ysollw(klon), ysolsw(klon)  
     REAL yfder(klon), ytaux(klon), ytauy(klon)  
     REAL yrugm(klon), yrads(klon), yrugoro(klon)  
   
166      REAL yfluxlat(klon)      REAL yfluxlat(klon)
   
167      REAL y_d_ts(klon)      REAL y_d_ts(klon)
168      REAL y_d_t(klon, klev), y_d_q(klon, klev)      REAL y_d_t(klon, klev), y_d_q(klon, klev)
169      REAL y_d_u(klon, klev), y_d_v(klon, klev)      REAL y_d_u(klon, klev), y_d_v(klon, klev)
170      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)      REAL y_flux_t(klon), y_flux_q(klon)
171      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)      REAL y_flux_u(klon), y_flux_v(klon)
172      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
173      REAL coefh(klon, klev), coefm(klon, klev)      REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
174        real ycdragh(klon), ycdragm(klon)
175      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
176      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
177      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)      REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
178        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)  
179      REAL delp(klon, klev)      REAL delp(klon, klev)
180      INTEGER i, k, nsrf      INTEGER i, k, nsrf
   
181      INTEGER ni(klon), knon, j      INTEGER ni(klon), knon, j
182    
183      REAL pctsrf_pot(klon, nbsrf)      REAL pctsrf_pot(klon, nbsrf)
184      ! "pourcentage potentiel" pour tenir compte des éventuelles      ! "pourcentage potentiel" pour tenir compte des \'eventuelles
185      ! apparitions ou disparitions de la glace de mer      ! apparitions ou disparitions de la glace de mer
186    
187      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.      REAL yt2m(klon), yq2m(klon), wind10m(klon)
188        REAL ustar(klon)
     ! maf pour sorties IOISPL en cas de debugagage  
   
     CHARACTER(80) cldebug  
     SAVE cldebug  
     CHARACTER(8) cl_surf(nbsrf)  
     SAVE cl_surf  
     INTEGER nhoridbg, nidbg  
     SAVE nhoridbg, nidbg  
     INTEGER ndexbg(iim*(jjm+1))  
     REAL zx_lon(iim, jjm+1), zx_lat(iim, jjm+1), zjulian  
     REAL tabindx(klon)  
     REAL debugtab(iim, jjm+1)  
     LOGICAL first_appel  
     SAVE first_appel  
     DATA first_appel/ .TRUE./  
     LOGICAL:: debugindex = .FALSE.  
     INTEGER idayref  
   
     REAL yt2m(klon), yq2m(klon), yu10m(klon)  
     REAL yustar(klon)  
     ! -- LOOP  
     REAL yu10mx(klon)  
     REAL yu10my(klon)  
     REAL ywindsp(klon)  
     ! -- LOOP  
189    
190      REAL yt10m(klon), yq10m(klon)      REAL yt10m(klon), yq10m(klon)
191      REAL ypblh(klon)      REAL ypblh(klon)
# Line 260  contains Line 195  contains
195      REAL ycteicl(klon)      REAL ycteicl(klon)
196      REAL ypblt(klon)      REAL ypblt(klon)
197      REAL ytherm(klon)      REAL ytherm(klon)
198      REAL ytrmb1(klon)      REAL u1(klon), v1(klon)
     REAL ytrmb2(klon)  
     REAL ytrmb3(klon)  
     REAL uzon(klon), vmer(klon)  
199      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
200      REAL psfce(klon), patm(klon)      REAL psfce(klon), patm(klon)
201        REAL zgeo1(klon)
     REAL qairsol(klon), zgeo1(klon)  
202      REAL rugo1(klon)      REAL rugo1(klon)
203        REAL zgeop(klon, klev)
     ! utiliser un jeu de fonctions simples                
     LOGICAL zxli  
     PARAMETER (zxli=.FALSE.)  
204    
205      !------------------------------------------------------------      !------------------------------------------------------------
206    
207      ytherm = 0.      albsol = sum(falbe * pctsrf, dim = 2)
208    
209      IF (debugindex .AND. first_appel) THEN      ! R\'epartition sous maille des flux longwave et shortwave
210         first_appel = .FALSE.      ! R\'epartition du longwave par sous-surface lin\'earis\'ee
211    
212         ! initialisation sorties netcdf      forall (nsrf = 1:nbsrf)
213           fsollw(:, nsrf) = sollw + 4. * RSIGMA * tsol**3 &
214                * (tsol - ftsol(:, nsrf))
215           fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
216        END forall
217    
218         idayref = day_ini      ytherm = 0.
        CALL ymds2ju(annee_ref, 1, idayref, 0., zjulian)  
        CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlon, zx_lon)  
        DO i = 1, iim  
           zx_lon(i, 1) = rlon(i+1)  
           zx_lon(i, jjm+1) = rlon(i+1)  
        END DO  
        CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlat, zx_lat)  
        cldebug = 'sous_index'  
        CALL histbeg_totreg(cldebug, zx_lon(:, 1), zx_lat(1, :), 1, &  
             iim, 1, jjm+1, itau_phy, zjulian, dtime, nhoridbg, nidbg)  
        ! no vertical axis  
        cl_surf(1) = 'ter'  
        cl_surf(2) = 'lic'  
        cl_surf(3) = 'oce'  
        cl_surf(4) = 'sic'  
        DO nsrf = 1, nbsrf  
           CALL histdef(nidbg, cl_surf(nsrf), cl_surf(nsrf), '-', iim, jjm+1, &  
                nhoridbg, 1, 1, 1, -99, 'inst', dtime, dtime)  
        END DO  
        CALL histend(nidbg)  
        CALL histsync(nidbg)  
     END IF  
219    
220      DO k = 1, klev ! epaisseur de couche      DO k = 1, klev ! epaisseur de couche
221         DO i = 1, klon         DO i = 1, klon
222            delp(i, k) = paprs(i, k) - paprs(i, k+1)            delp(i, k) = paprs(i, k) - paprs(i, k + 1)
223         END DO         END DO
224      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  
225    
226      ! Initialization:      ! Initialization:
227      rugmer = 0.      rugmer = 0.
# Line 325  contains Line 229  contains
229      cdragm = 0.      cdragm = 0.
230      dflux_t = 0.      dflux_t = 0.
231      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.  
     ytaux = 0.  
     ytauy = 0.  
     ysolsw = 0.  
     ysollw = 0.  
232      yrugos = 0.      yrugos = 0.
     yu1 = 0.  
     yv1 = 0.  
     yrads = 0.  
233      ypaprs = 0.      ypaprs = 0.
234      ypplay = 0.      ypplay = 0.
235      ydelp = 0.      ydelp = 0.
     yu = 0.  
     yv = 0.  
     yt = 0.  
     yq = 0.  
     pctsrf_new = 0.  
     y_flux_u = 0.  
     y_flux_v = 0.  
     !$$ PB  
     y_dflux_t = 0.  
     y_dflux_q = 0.  
     ytsoil = 999999.  
236      yrugoro = 0.      yrugoro = 0.
     ! -- LOOP  
     yu10mx = 0.  
     yu10my = 0.  
     ywindsp = 0.  
     ! -- LOOP  
237      d_ts = 0.      d_ts = 0.
     !§§§ PB  
     yfluxlat = 0.  
238      flux_t = 0.      flux_t = 0.
239      flux_q = 0.      flux_q = 0.
240      flux_u = 0.      flux_u = 0.
241      flux_v = 0.      flux_v = 0.
242        fluxlat = 0.
243      d_t = 0.      d_t = 0.
244      d_q = 0.      d_q = 0.
245      d_u = 0.      d_u = 0.
246      d_v = 0.      d_v = 0.
247      ycoefh = 0.      coefh = 0.
248        fqcalving = 0.
249        run_off_lic = 0.
250    
251        ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
252        ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
253        ! (\`a affiner).
254    
255        pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
256        pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
257        pctsrf_pot(:, is_oce) = 1. - masque
258        pctsrf_pot(:, is_sic) = 1. - masque
259    
260        ! Tester si c'est le moment de lire le fichier:
261        if (mod(itap - 1, lmt_pas) == 0) then
262           CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
263        endif
264    
265      ! Boucler sur toutes les sous-fractions du sol:      ! Boucler sur toutes les sous-fractions du sol:
266    
     ! Initialisation des "pourcentages potentiels". On considère ici qu'on  
     ! peut avoir potentiellement de la glace sur tout le domaine océanique  
     ! (à affiner)  
   
     pctsrf_pot = pctsrf  
     pctsrf_pot(:, is_oce) = 1. - zmasq  
     pctsrf_pot(:, is_sic) = 1. - zmasq  
   
267      loop_surface: DO nsrf = 1, nbsrf      loop_surface: DO nsrf = 1, nbsrf
268         ! Chercher les indices :         ! Define ni and knon:
269    
270         ni = 0         ni = 0
271         knon = 0         knon = 0
272    
273         DO i = 1, klon         DO i = 1, klon
274            ! Pour déterminer le domaine à traiter, on utilise les surfaces            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
275            ! "potentielles"            ! "potentielles"
276            IF (pctsrf_pot(i, nsrf) > epsfra) THEN            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
277               knon = knon + 1               knon = knon + 1
# Line 400  contains Line 279  contains
279            END IF            END IF
280         END DO         END DO
281    
        ! variables pour avoir une sortie IOIPSL des INDEX  
        IF (debugindex) THEN  
           tabindx = 0.  
           DO i = 1, knon  
              tabindx(i) = real(i)  
           END DO  
           debugtab = 0.  
           ndexbg = 0  
           CALL gath2cpl(tabindx, debugtab, klon, knon, iim, jjm, ni)  
           CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)  
        END IF  
   
282         if_knon: IF (knon /= 0) then         if_knon: IF (knon /= 0) then
283            DO j = 1, knon            ypctsrf(:knon) = pctsrf(ni(:knon), nsrf)
284               i = ni(j)            yts(:knon) = ftsol(ni(:knon), nsrf)
285               ypct(j) = pctsrf(i, nsrf)            snow(:knon) = fsnow(ni(:knon), nsrf)
286               yts(j) = ts(i, nsrf)            yqsurf(:knon) = fqsurf(ni(:knon), nsrf)
287               ytslab(i) = tslab(i)            yalb(:knon) = falbe(ni(:knon), nsrf)
288               ysnow(j) = snow(i, nsrf)            yrain_fall(:knon) = rain_fall(ni(:knon))
289               yqsurf(j) = qsurf(i, nsrf)            ysnow_fall(:knon) = snow_fall(ni(:knon))
290               yalb(j) = albe(i, nsrf)            yagesno(:knon) = agesno(ni(:knon), nsrf)
291               yalblw(j) = alblw(i, nsrf)            yrugos(:knon) = frugs(ni(:knon), nsrf)
292               yrain_f(j) = rain_fall(i)            yrugoro(:knon) = rugoro(ni(:knon))
293               ysnow_f(j) = snow_f(i)            radsol(:knon) = fsolsw(ni(:knon), nsrf) + fsollw(ni(:knon), nsrf)
294               yagesno(j) = agesno(i, nsrf)            ypaprs(:knon, klev + 1) = paprs(ni(:knon), klev + 1)
295               yfder(j) = fder(i)            y_run_off_lic_0(:knon) = run_off_lic_0(ni(:knon))
              ytaux(j) = flux_u(i, 1, nsrf)  
              ytauy(j) = flux_v(i, 1, nsrf)  
              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  
296    
297            ! IF bucket model for continent, copy soil water content            ! For continent, copy soil water content
298            IF (nsrf == is_ter .AND. .NOT. ok_veget) THEN            IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
              DO j = 1, knon  
                 i = ni(j)  
                 yqsol(j) = qsol(i)  
              END DO  
           ELSE  
              yqsol = 0.  
           END IF  
299    
300            DO k = 1, nsoilmx            ytsoil(:knon, :) = ftsoil(ni(:knon), :, nsrf)
              DO j = 1, knon  
                 i = ni(j)  
                 ytsoil(j, k) = ftsoil(i, k, nsrf)  
              END DO  
           END DO  
301    
302            DO k = 1, klev            DO k = 1, klev
303               DO j = 1, knon               DO j = 1, knon
304                  i = ni(j)                  i = ni(j)
305                  ypaprs(j, k) = paprs(i, k)                  ypaprs(j, k) = paprs(i, k)
306                  ypplay(j, k) = pplay(i, k)                  ypplay(j, k) = play(i, k)
307                  ydelp(j, k) = delp(i, k)                  ydelp(j, k) = delp(i, k)
308                  yu(j, k) = u(i, k)                  yu(j, k) = u(i, k)
309                  yv(j, k) = v(i, k)                  yv(j, k) = v(i, k)
# Line 472  contains Line 312  contains
312               END DO               END DO
313            END DO            END DO
314    
315            ! calculer Cdrag et les coefficients d'echange            ! Calculer les géopotentiels de chaque couche:
           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  
316    
317            ! on met un seuil pour coefm et coefh            zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
318            IF (nsrf == is_oce) THEN                 + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
              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  
   
           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  
              DO k = 1, klev + 1  
                 DO j = 1, knon  
                    i = ni(j)  
                    yq2(j, k) = q2(i, k, nsrf)  
                 END DO  
              END DO  
   
              CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)  
   
              IF (prt_level > 9) THEN  
                 PRINT *, 'USTAR = ', yustar  
              END IF  
319    
320               ! iflag_pbl peut être utilisé comme longueur de mélange            DO k = 2, klev
321                 zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
322                      * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
323                      * (ypplay(:knon, k - 1) - ypplay(:knon, k))
324              ENDDO
325    
326              CALL cdrag(nsrf, sqrt(yu(:knon, 1)**2 + yv(:knon, 1)**2), &
327                   yt(:knon, 1), yq(:knon, 1), zgeop(:knon, 1), ypaprs(:knon, 1), &
328                   yts(:knon), yqsurf(:knon), yrugos(:knon), ycdragm(:knon), &
329                   ycdragh(:knon))
330    
331               IF (iflag_pbl >= 11) THEN            IF (iflag_pbl == 1) THEN
332                  CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &               ycdragm(:knon) = max(ycdragm(:knon), 0.)
333                       yu, yv, yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, &               ycdragh(:knon) = max(ycdragh(:knon), 0.)
334                       yustar, 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  
335    
336               coefm(:knon, 2:) = ykmm(:knon, 2:klev)            ! on met un seuil pour ycdragm et ycdragh
337               coefh(:knon, 2:) = ykmn(:knon, 2:klev)            IF (nsrf == is_oce) THEN
338                 ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
339                 ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
340            END IF            END IF
341    
342            ! calculer la diffusion des vitesses "u" et "v"            IF (iflag_pbl >= 6) yq2(:knon, :) = q2(ni(:knon), :, nsrf)
343            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &            call coef_diff_turb(nsrf, ni(:knon), ypaprs(:knon, :), &
344                 ypplay, ydelp, y_d_u, y_flux_u)                 ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
345            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &                 yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
346                 ypplay, ydelp, y_d_v, y_flux_v)                 ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
347    
348            ! pour le couplage            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
349            ytaux = y_flux_u(:, 1)                 ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
350            ytauy = y_flux_v(:, 1)                 ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
351                   y_flux_u(:knon))
352            ! calculer la diffusion de "q" et de "h"            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
353            CALL clqh(dtime, itap, jour, debut, rlat, knon, nsrf, ni, pctsrf, &                 ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
354                 soil_model, ytsoil, yqsol, ok_veget, ocean, rmu0, co2_ppm, &                 ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
355                 yrugos, yrugoro, yu1, yv1, coefh(:knon, :), yt, yq, yts, &                 y_flux_v(:knon))
356                 ypaprs, ypplay, ydelp, yrads, yalb, yalblw, ysnow, yqsurf, &  
357                 yrain_f, ysnow_f, yfder, ysolsw, yfluxlat, pctsrf_new, &            CALL clqh(julien, nsrf, ni(:knon), ytsoil(:knon, :), yqsol(:knon), &
358                 yagesno, y_d_t, y_d_q, y_d_ts, yz0_new, y_flux_t, y_flux_q, &                 mu0(ni(:knon)), yrugos(:knon), yrugoro(:knon), yu(:knon, 1), &
359                 y_dflux_t, y_dflux_q, y_fqcalving, y_ffonte, y_run_off_lic_0, &                 yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), yt(:knon, :), &
360                 y_flux_o, y_flux_g, ytslab, y_seaice)                 yq(:knon, :), yts(:knon), ypaprs(:knon, :), ypplay(:knon, :), &
361                   ydelp(:knon, :), radsol(:knon), yalb(:knon), snow(:knon), &
362                   yqsurf(:knon), yrain_fall(:knon), ysnow_fall(:knon), &
363                   yfluxlat(:knon), pctsrf_new_sic(ni(:knon)), yagesno(:knon), &
364                   y_d_t(:knon, :), y_d_q(:knon, :), y_d_ts(:knon), &
365                   yz0_new(:knon), y_flux_t(:knon), y_flux_q(:knon), &
366                   y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving(:knon), &
367                   y_ffonte(:knon), y_run_off_lic_0(:knon), y_run_off_lic(:knon))
368    
369            ! calculer la longueur de rugosite sur ocean            ! calculer la longueur de rugosite sur ocean
370    
371            yrugm = 0.            yrugm = 0.
372    
373            IF (nsrf == is_oce) THEN            IF (nsrf == is_oce) THEN
374               DO j = 1, knon               DO j = 1, knon
375                  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) &
376                       0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))                       / rg + 0.11 * 14E-6 &
377                         / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
378                  yrugm(j) = max(1.5E-05, yrugm(j))                  yrugm(j) = max(1.5E-05, yrugm(j))
379               END DO               END DO
380            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  
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) * ypctsrf(j)
386                  coefm(j, k) = coefm(j, k)*ypct(j)                  y_d_q(j, k) = y_d_q(j, k) * ypctsrf(j)
387                  y_d_t(j, k) = y_d_t(j, k)*ypct(j)                  y_d_u(j, k) = y_d_u(j, k) * ypctsrf(j)
388                  y_d_q(j, k) = y_d_q(j, k)*ypct(j)                  y_d_v(j, k) = y_d_v(j, k) * ypctsrf(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               DO j = 1, knon               qsol(ni(:knon)) = yqsol(:knon)
423                  i = ni(j)            else IF (nsrf == is_lic) THEN
                 qsol(i) = yqsol(j)  
              END DO  
           END IF  
           IF (nsrf == is_lic) THEN  
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            !$$$ PB ajout pour soil  
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 656  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            !cc 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              ! diagnostic t, q a 2m et u, v a 10m
448    
449            DO j = 1, knon            DO j = 1, knon
450               i = ni(j)               i = ni(j)
451               uzon(j) = yu(j, 1) + y_d_u(j, 1)               u1(j) = yu(j, 1) + y_d_u(j, 1)
452               vmer(j) = yv(j, 1) + y_d_v(j, 1)               v1(j) = yv(j, 1) + y_d_v(j, 1)
453               tair1(j) = yt(j, 1) + y_d_t(j, 1)               tair1(j) = yt(j, 1) + y_d_t(j, 1)
454               qair1(j) = yq(j, 1) + y_d_q(j, 1)               qair1(j) = yq(j, 1) + y_d_q(j, 1)
455               zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &               zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
456                    1)))*(ypaprs(j, 1)-ypplay(j, 1))                    1))) * (ypaprs(j, 1)-ypplay(j, 1))
457               tairsol(j) = yts(j) + y_d_ts(j)               tairsol(j) = yts(j) + y_d_ts(j)
458               rugo1(j) = yrugos(j)               rugo1(j) = yrugos(j)
459               IF (nsrf == is_oce) THEN               IF (nsrf == is_oce) THEN
460                  rugo1(j) = rugos(i, nsrf)                  rugo1(j) = frugs(i, nsrf)
461               END IF               END IF
462               psfce(j) = ypaprs(j, 1)               psfce(j) = ypaprs(j, 1)
463               patm(j) = ypplay(j, 1)               patm(j) = ypplay(j, 1)
   
              qairsol(j) = yqsurf(j)  
464            END DO            END DO
465    
466            CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &            CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), qair1, &
467                 zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &                 zgeo1, tairsol, yqsurf(:knon), rugo1, psfce, patm, yt2m, yq2m, &
468                 yt10m, yq10m, yu10m, yustar)                 yt10m, yq10m, wind10m(:knon), ustar(:knon))
469    
470            DO j = 1, knon            DO j = 1, knon
471               i = ni(j)               i = ni(j)
472               t2m(i, nsrf) = yt2m(j)               t2m(i, nsrf) = yt2m(j)
473               q2m(i, nsrf) = yq2m(j)               q2m(i, nsrf) = yq2m(j)
474    
475               ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman               u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
476               u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)                    / sqrt(u1(j)**2 + v1(j)**2)
477               v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)               v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
478                      / sqrt(u1(j)**2 + v1(j)**2)
479            END DO            END DO
480    
481            CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
482                 y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &                 y_flux_q(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
483                 ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)                 yq(:knon, :), ypblh(:knon), ycapcl, yoliqcl, ycteicl, ypblt, &
484                   ytherm, ylcl)
485    
486            DO j = 1, knon            DO j = 1, knon
487               i = ni(j)               i = ni(j)
# Line 709  contains Line 492  contains
492               cteicl(i, nsrf) = ycteicl(j)               cteicl(i, nsrf) = ycteicl(j)
493               pblt(i, nsrf) = ypblt(j)               pblt(i, nsrf) = ypblt(j)
494               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)  
495            END DO            END DO
496    
497            DO j = 1, knon            IF (iflag_pbl >= 6) q2(ni(:knon), :, nsrf) = yq2(:knon, :)
498               DO k = 1, klev + 1         else
499                  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  
           IF (ocean == 'slab  ') THEN  
              IF (nsrf == is_oce) THEN  
                 tslab(1:klon) = ytslab(1:klon)  
                 seaice(1:klon) = y_seaice(1:klon)  
              END IF  
           END IF  
500         end IF if_knon         end IF if_knon
501      END DO loop_surface      END DO loop_surface
502    
503      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
504        frugs(:, is_oce) = rugmer
505        pctsrf(:, is_oce) = pctsrf_new_oce
506        pctsrf(:, is_sic) = pctsrf_new_sic
507    
508        CALL histwrite_phy("run_off_lic", run_off_lic)
509        ftsol = ftsol + d_ts ! update surface temperature
510        CALL histwrite_phy("dtsvdfo", d_ts(:, is_oce))
511        CALL histwrite_phy("dtsvdft", d_ts(:, is_ter))
512        CALL histwrite_phy("dtsvdfg", d_ts(:, is_lic))
513        CALL histwrite_phy("dtsvdfi", d_ts(:, is_sic))
514    
515      rugos(:, is_oce) = rugmer    END SUBROUTINE pbl_surface
     pctsrf = pctsrf_new  
   
   END SUBROUTINE clmain  
516    
517  end module clmain_m  end module pbl_surface_m
Legend:
Removed from v.76  
changed lines
  Added in v.328
  ViewVC Help
Powered by ViewVC 1.1.21