/[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/libf/phylmd/clmain.f90 revision 72 by guez, Tue Jul 23 13:00:07 2013 UTC trunk/phylmd/Interface_surf/pbl_surface.f revision 310 by guez, Thu Sep 27 16:29:06 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, &    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, d_ts, &
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: zmasq
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
41    
     ! Arguments:  
   
     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        ! 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 co2_ppm ! taux CO2 atmosphere  
51      LOGICAL ok_veget      REAL, INTENT(IN):: 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      LOGICAL, INTENT(IN):: debut      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        REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)
68    
69        REAL, intent(in):: rain_fall(klon)
70        ! liquid water mass flux (kg / m2 / s), positive down
71    
72        REAL, intent(in):: snow_fall(klon)
73        ! solid water mass flux (kg / m2 / s), positive down
74    
75        REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)
76      real agesno(klon, nbsrf)      real agesno(klon, nbsrf)
77      REAL, INTENT(IN):: rugoro(klon)      REAL, INTENT(IN):: rugoro(klon)
78    
79      REAL d_t(klon, klev), d_q(klon, klev)      REAL, intent(out):: d_t(:, :), d_q(:, :) ! (klon, klev)
80      ! d_t------output-R- le changement pour "t"      ! changement pour t et q
     ! d_q------output-R- le changement pour "q"  
81    
82      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
83      ! changement pour "u" et "v"      ! changement pour "u" et "v"
84    
85      REAL d_ts(klon, nbsrf)      REAL, intent(out):: d_ts(:, :) ! (klon, nbsrf) variation of ftsol
     ! d_ts-----output-R- le changement pour "ts"  
86    
87      REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)      REAL, intent(out):: flux_t(klon, nbsrf)
88      ! 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
89      !                    (orientation positive vers le bas)      ! vers le bas) à la surface
90      ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)  
91        REAL, intent(out):: flux_q(klon, nbsrf)
92      REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)      ! flux de vapeur d'eau (kg / m2 / s) à la surface
93      ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal  
94      ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal      REAL, intent(out):: flux_u(:, :), flux_v(:, :) ! (klon, nbsrf)
95        ! tension du vent (flux turbulent de vent) à la surface, en Pa
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 fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
146      real ytslab(klon)      REAL fsolsw(klon, nbsrf) ! flux solaire absorb\'e pour chaque sous-surface
     real y_seaice(klon)  
     REAL y_fqcalving(klon), y_ffonte(klon)  
     real y_run_off_lic_0(klon)  
147    
148      REAL rugmer(klon)      ! la nouvelle repartition des surfaces sortie de l'interface
149        REAL, save:: pctsrf_new_oce(klon)
150        REAL, save:: pctsrf_new_sic(klon)
151    
152        REAL y_fqcalving(klon), y_ffonte(klon)
153        real y_run_off_lic_0(klon), y_run_off_lic(klon)
154        REAL run_off_lic(klon) ! ruissellement total
155        REAL rugmer(klon)
156      REAL ytsoil(klon, nsoilmx)      REAL ytsoil(klon, nsoilmx)
157        REAL yts(klon), ypctsrf(klon), yz0_new(klon)
158      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      real yrugos(klon) ! longueur de rugosite (en m)
159      REAL yalb(klon)      REAL yalb(klon)
160      REAL yalblw(klon)      REAL snow(klon), yqsurf(klon), yagesno(klon)
161      REAL yu1(klon), yv1(klon)      real yqsol(klon) ! column-density of water in soil, in kg m-2
162      ! 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
163      ! la premiere couche      REAL ysnow_fall(klon) ! solid water mass flux (kg / m2 / s), positive down
164      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)  
   
165      REAL yfluxlat(klon)      REAL yfluxlat(klon)
   
166      REAL y_d_ts(klon)      REAL y_d_ts(klon)
167      REAL y_d_t(klon, klev), y_d_q(klon, klev)      REAL y_d_t(klon, klev), y_d_q(klon, klev)
168      REAL y_d_u(klon, klev), y_d_v(klon, klev)      REAL y_d_u(klon, klev), y_d_v(klon, klev)
169      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)      REAL y_flux_t(klon), y_flux_q(klon)
170      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)      REAL y_flux_u(klon), y_flux_v(klon)
171      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
172      REAL coefh(klon, klev), coefm(klon, klev)      REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
173        real ycdragh(klon), ycdragm(klon)
174      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
175      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
176      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)      REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
177        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)  
178      REAL delp(klon, klev)      REAL delp(klon, klev)
179      INTEGER i, k, nsrf      INTEGER i, k, nsrf
   
180      INTEGER ni(klon), knon, j      INTEGER ni(klon), knon, j
181    
182      REAL pctsrf_pot(klon, nbsrf)      REAL pctsrf_pot(klon, nbsrf)
183      ! "pourcentage potentiel" pour tenir compte des éventuelles      ! "pourcentage potentiel" pour tenir compte des \'eventuelles
184      ! apparitions ou disparitions de la glace de mer      ! apparitions ou disparitions de la glace de mer
185    
186      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.      REAL yt2m(klon), yq2m(klon), wind10m(klon)
187        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  
188    
189      REAL yt10m(klon), yq10m(klon)      REAL yt10m(klon), yq10m(klon)
190      REAL ypblh(klon)      REAL ypblh(klon)
# Line 260  contains Line 194  contains
194      REAL ycteicl(klon)      REAL ycteicl(klon)
195      REAL ypblt(klon)      REAL ypblt(klon)
196      REAL ytherm(klon)      REAL ytherm(klon)
197      REAL ytrmb1(klon)      REAL u1(klon), v1(klon)
     REAL ytrmb2(klon)  
     REAL ytrmb3(klon)  
     REAL uzon(klon), vmer(klon)  
198      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
199      REAL psfce(klon), patm(klon)      REAL psfce(klon), patm(klon)
200        REAL zgeo1(klon)
     REAL qairsol(klon), zgeo1(klon)  
201      REAL rugo1(klon)      REAL rugo1(klon)
202        REAL zgeop(klon, klev)
     ! utiliser un jeu de fonctions simples                
     LOGICAL zxli  
     PARAMETER (zxli=.FALSE.)  
203    
204      !------------------------------------------------------------      !------------------------------------------------------------
205    
206      ytherm = 0.      albsol = sum(falbe * pctsrf, dim = 2)
207    
208      IF (debugindex .AND. first_appel) THEN      ! R\'epartition sous maille des flux longwave et shortwave
209         first_appel = .FALSE.      ! R\'epartition du longwave par sous-surface lin\'earis\'ee
210    
211         ! initialisation sorties netcdf      forall (nsrf = 1:nbsrf)
212           fsollw(:, nsrf) = sollw + 4. * RSIGMA * tsol**3 &
213                * (tsol - ftsol(:, nsrf))
214           fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
215        END forall
216    
217         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  
218    
219      DO k = 1, klev ! epaisseur de couche      DO k = 1, klev ! epaisseur de couche
220         DO i = 1, klon         DO i = 1, klon
221            delp(i, k) = paprs(i, k) - paprs(i, k+1)            delp(i, k) = paprs(i, k) - paprs(i, k + 1)
222         END DO         END DO
223      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  
224    
225      ! Initialization:      ! Initialization:
226      rugmer = 0.      rugmer = 0.
# Line 325  contains Line 228  contains
228      cdragm = 0.      cdragm = 0.
229      dflux_t = 0.      dflux_t = 0.
230      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.  
231      yrugos = 0.      yrugos = 0.
     yu1 = 0.  
     yv1 = 0.  
     yrads = 0.  
232      ypaprs = 0.      ypaprs = 0.
233      ypplay = 0.      ypplay = 0.
234      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.  
235      yrugoro = 0.      yrugoro = 0.
     ! -- LOOP  
     yu10mx = 0.  
     yu10my = 0.  
     ywindsp = 0.  
     ! -- LOOP  
236      d_ts = 0.      d_ts = 0.
     !§§§ PB  
     yfluxlat = 0.  
237      flux_t = 0.      flux_t = 0.
238      flux_q = 0.      flux_q = 0.
239      flux_u = 0.      flux_u = 0.
240      flux_v = 0.      flux_v = 0.
241        fluxlat = 0.
242      d_t = 0.      d_t = 0.
243      d_q = 0.      d_q = 0.
244      d_u = 0.      d_u = 0.
245      d_v = 0.      d_v = 0.
246      ycoefh = 0.      coefh = 0.
247        fqcalving = 0.
248      ! Boucler sur toutes les sous-fractions du sol:      run_off_lic = 0.
249    
250      ! Initialisation des "pourcentages potentiels". On considère ici qu'on      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
251      ! peut avoir potentiellement de la glace sur tout le domaine océanique      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
252      ! (à affiner)      ! (\`a affiner).
253    
254      pctsrf_pot = pctsrf      pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
255        pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
256      pctsrf_pot(:, is_oce) = 1. - zmasq      pctsrf_pot(:, is_oce) = 1. - zmasq
257      pctsrf_pot(:, is_sic) = 1. - zmasq      pctsrf_pot(:, is_sic) = 1. - zmasq
258    
259        ! Tester si c'est le moment de lire le fichier:
260        if (mod(itap - 1, lmt_pas) == 0) then
261           CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
262        endif
263    
264        ! Boucler sur toutes les sous-fractions du sol:
265    
266      loop_surface: DO nsrf = 1, nbsrf      loop_surface: DO nsrf = 1, nbsrf
267         ! Chercher les indices :         ! Define ni and knon:
268    
269         ni = 0         ni = 0
270         knon = 0         knon = 0
271    
272         DO i = 1, klon         DO i = 1, klon
273            ! Pour déterminer le domaine à traiter, on utilise les surfaces            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
274            ! "potentielles"            ! "potentielles"
275            IF (pctsrf_pot(i, nsrf) > epsfra) THEN            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
276               knon = knon + 1               knon = knon + 1
# Line 400  contains Line 278  contains
278            END IF            END IF
279         END DO         END DO
280    
        ! 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  
   
281         if_knon: IF (knon /= 0) then         if_knon: IF (knon /= 0) then
282            DO j = 1, knon            ypctsrf(:knon) = pctsrf(ni(:knon), nsrf)
283               i = ni(j)            yts(:knon) = ftsol(ni(:knon), nsrf)
284               ypct(j) = pctsrf(i, nsrf)            snow(:knon) = fsnow(ni(:knon), nsrf)
285               yts(j) = ts(i, nsrf)            yqsurf(:knon) = fqsurf(ni(:knon), nsrf)
286               ytslab(i) = tslab(i)            yalb(:knon) = falbe(ni(:knon), nsrf)
287               ysnow(j) = snow(i, nsrf)            yrain_fall(:knon) = rain_fall(ni(:knon))
288               yqsurf(j) = qsurf(i, nsrf)            ysnow_fall(:knon) = snow_fall(ni(:knon))
289               yalb(j) = albe(i, nsrf)            yagesno(:knon) = agesno(ni(:knon), nsrf)
290               yalblw(j) = alblw(i, nsrf)            yrugos(:knon) = frugs(ni(:knon), nsrf)
291               yrain_f(j) = rain_fall(i)            yrugoro(:knon) = rugoro(ni(:knon))
292               ysnow_f(j) = snow_f(i)            radsol(:knon) = fsolsw(ni(:knon), nsrf) + fsollw(ni(:knon), nsrf)
293               yagesno(j) = agesno(i, nsrf)            ypaprs(:knon, klev + 1) = paprs(ni(:knon), klev + 1)
294               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  
295    
296            ! IF bucket model for continent, copy soil water content            ! For continent, copy soil water content
297            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  
298    
299            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  
300    
301            DO k = 1, klev            DO k = 1, klev
302               DO j = 1, knon               DO j = 1, knon
303                  i = ni(j)                  i = ni(j)
304                  ypaprs(j, k) = paprs(i, k)                  ypaprs(j, k) = paprs(i, k)
305                  ypplay(j, k) = pplay(i, k)                  ypplay(j, k) = play(i, k)
306                  ydelp(j, k) = delp(i, k)                  ydelp(j, k) = delp(i, k)
307                  yu(j, k) = u(i, k)                  yu(j, k) = u(i, k)
308                  yv(j, k) = v(i, k)                  yv(j, k) = v(i, k)
# Line 472  contains Line 311  contains
311               END DO               END DO
312            END DO            END DO
313    
314            ! calculer Cdrag et les coefficients d'echange            ! Calculer les géopotentiels de chaque couche:
           CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &  
                yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))  
           IF (iflag_pbl == 1) THEN  
              CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)  
              coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))  
              coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))  
           END IF  
   
           ! on met un seuil pour coefm et coefh  
           IF (nsrf == is_oce) THEN  
              coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)  
              coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)  
           END IF  
   
           IF (ok_kzmin) THEN  
              ! Calcul d'une diffusion minimale pour les conditions tres stables  
              CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &  
                   coefm(:knon, 1), ycoefm0, ycoefh0)  
              coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))  
              coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))  
            END IF  
   
           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)  
315    
316               IF (prt_level > 9) THEN            zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
317                  PRINT *, 'USTAR = ', yustar                 + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
              END IF  
318    
319               ! iflag_pbl peut être utilisé comme longueur de mélange            DO k = 2, klev
320                 zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
321                      * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
322                      * (ypplay(:knon, k - 1) - ypplay(:knon, k))
323              ENDDO
324    
325              CALL cdrag(nsrf, sqrt(yu(:knon, 1)**2 + yv(:knon, 1)**2), &
326                   yt(:knon, 1), yq(:knon, 1), zgeop(:knon, 1), ypaprs(:knon, 1), &
327                   yts(:knon), yqsurf(:knon), yrugos(:knon), ycdragm(:knon), &
328                   ycdragh(:knon))
329    
330               IF (iflag_pbl >= 11) THEN            IF (iflag_pbl == 1) THEN
331                  CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &               ycdragm(:knon) = max(ycdragm(:knon), 0.)
332                       yu, yv, yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, &               ycdragh(:knon) = max(ycdragh(:knon), 0.)
333                       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  
334    
335               coefm(:knon, 2:) = ykmm(:knon, 2:klev)            ! on met un seuil pour ycdragm et ycdragh
336               coefh(:knon, 2:) = ykmn(:knon, 2:klev)            IF (nsrf == is_oce) THEN
337                 ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
338                 ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
339            END IF            END IF
340    
341            ! calculer la diffusion des vitesses "u" et "v"            IF (iflag_pbl >= 6) yq2(:knon, :) = q2(ni(:knon), :, nsrf)
342            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &            call coef_diff_turb(nsrf, ni(:knon), ypaprs(:knon, :), &
343                 ypplay, ydelp, y_d_u, y_flux_u)                 ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
344            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &                 yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
345                 ypplay, ydelp, y_d_v, y_flux_v)                 ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
346    
347            ! pour le couplage            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
348            ytaux = y_flux_u(:, 1)                 ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
349            ytauy = y_flux_v(:, 1)                 ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
350                   y_flux_u(:knon))
351            ! calculer la diffusion de "q" et de "h"            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
352            CALL clqh(dtime, itap, jour, debut, rlat, knon, nsrf, ni, pctsrf, &                 ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
353                 soil_model, ytsoil, yqsol, ok_veget, ocean, rmu0, co2_ppm, &                 ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
354                 yrugos, yrugoro, yu1, yv1, coefh(:knon, :), yt, yq, yts, &                 y_flux_v(:knon))
355                 ypaprs, ypplay, ydelp, yrads, yalb, yalblw, ysnow, yqsurf, &  
356                 yrain_f, ysnow_f, yfder, ysolsw, yfluxlat, pctsrf_new, &            CALL clqh(julien, nsrf, ni(:knon), ytsoil(:knon, :), yqsol(:knon), &
357                 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), &
358                 y_dflux_t, y_dflux_q, y_fqcalving, y_ffonte, y_run_off_lic_0, &                 yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), yt(:knon, :), &
359                 y_flux_o, y_flux_g, ytslab, y_seaice)                 yq(:knon, :), yts(:knon), ypaprs(:knon, :), ypplay(:knon, :), &
360                   ydelp(:knon, :), radsol(:knon), yalb(:knon), snow(:knon), &
361                   yqsurf(:knon), yrain_fall(:knon), ysnow_fall(:knon), &
362                   yfluxlat(:knon), pctsrf_new_sic(ni(:knon)), yagesno(:knon), &
363                   y_d_t(:knon, :), y_d_q(:knon, :), y_d_ts(:knon), &
364                   yz0_new(:knon), y_flux_t(:knon), y_flux_q(:knon), &
365                   y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving(:knon), &
366                   y_ffonte(:knon), y_run_off_lic_0(:knon), y_run_off_lic(:knon))
367    
368            ! calculer la longueur de rugosite sur ocean            ! calculer la longueur de rugosite sur ocean
369    
370            yrugm = 0.            yrugm = 0.
371    
372            IF (nsrf == is_oce) THEN            IF (nsrf == is_oce) THEN
373               DO j = 1, knon               DO j = 1, knon
374                  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) &
375                       0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))                       / rg + 0.11 * 14E-6 &
376                         / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
377                  yrugm(j) = max(1.5E-05, yrugm(j))                  yrugm(j) = max(1.5E-05, yrugm(j))
378               END DO               END DO
379            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  
380    
381            DO k = 1, klev            DO k = 1, klev
382               DO j = 1, knon               DO j = 1, knon
383                  i = ni(j)                  i = ni(j)
384                  coefh(j, k) = coefh(j, k)*ypct(j)                  y_d_t(j, k) = y_d_t(j, k) * ypctsrf(j)
385                  coefm(j, k) = coefm(j, k)*ypct(j)                  y_d_q(j, k) = y_d_q(j, k) * ypctsrf(j)
386                  y_d_t(j, k) = y_d_t(j, k)*ypct(j)                  y_d_u(j, k) = y_d_u(j, k) * ypctsrf(j)
387                  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)  
388               END DO               END DO
389            END DO            END DO
390    
391            evap(:, nsrf) = -flux_q(:, 1, nsrf)            flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
392              flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
393            albe(:, nsrf) = 0.            flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
394            alblw(:, nsrf) = 0.            flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
395            snow(:, nsrf) = 0.  
396            qsurf(:, nsrf) = 0.            falbe(:, nsrf) = 0.
397            rugos(:, nsrf) = 0.            fsnow(:, nsrf) = 0.
398            fluxlat(:, nsrf) = 0.            fqsurf(:, nsrf) = 0.
399              frugs(:, nsrf) = 0.
400            DO j = 1, knon            DO j = 1, knon
401               i = ni(j)               i = ni(j)
402               d_ts(i, nsrf) = y_d_ts(j)               d_ts(i, nsrf) = y_d_ts(j)
403               albe(i, nsrf) = yalb(j)               falbe(i, nsrf) = yalb(j)
404               alblw(i, nsrf) = yalblw(j)               fsnow(i, nsrf) = snow(j)
405               snow(i, nsrf) = ysnow(j)               fqsurf(i, nsrf) = yqsurf(j)
406               qsurf(i, nsrf) = yqsurf(j)               frugs(i, nsrf) = yz0_new(j)
              rugos(i, nsrf) = yz0_new(j)  
407               fluxlat(i, nsrf) = yfluxlat(j)               fluxlat(i, nsrf) = yfluxlat(j)
408               IF (nsrf == is_oce) THEN               IF (nsrf == is_oce) THEN
409                  rugmer(i) = yrugm(j)                  rugmer(i) = yrugm(j)
410                  rugos(i, nsrf) = yrugm(j)                  frugs(i, nsrf) = yrugm(j)
411               END IF               END IF
412               agesno(i, nsrf) = yagesno(j)               agesno(i, nsrf) = yagesno(j)
413               fqcalving(i, nsrf) = y_fqcalving(j)               fqcalving(i, nsrf) = y_fqcalving(j)
414               ffonte(i, nsrf) = y_ffonte(j)               ffonte(i, nsrf) = y_ffonte(j)
415               cdragh(i) = cdragh(i) + coefh(j, 1)               cdragh(i) = cdragh(i) + ycdragh(j) * ypctsrf(j)
416               cdragm(i) = cdragm(i) + coefm(j, 1)               cdragm(i) = cdragm(i) + ycdragm(j) * ypctsrf(j)
417               dflux_t(i) = dflux_t(i) + y_dflux_t(j)               dflux_t(i) = dflux_t(i) + y_dflux_t(j) * ypctsrf(j)
418               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)  
419            END DO            END DO
420            IF (nsrf == is_ter) THEN            IF (nsrf == is_ter) THEN
421               DO j = 1, knon               qsol(ni(:knon)) = yqsol(:knon)
422                  i = ni(j)            else IF (nsrf == is_lic) THEN
                 qsol(i) = yqsol(j)  
              END DO  
           END IF  
           IF (nsrf == is_lic) THEN  
423               DO j = 1, knon               DO j = 1, knon
424                  i = ni(j)                  i = ni(j)
425                  run_off_lic_0(i) = y_run_off_lic_0(j)                  run_off_lic_0(i) = y_run_off_lic_0(j)
426                    run_off_lic(i) = y_run_off_lic(j)
427               END DO               END DO
428            END IF            END IF
429            !$$$ PB ajout pour soil  
430            ftsoil(:, :, nsrf) = 0.            ftsoil(:, :, nsrf) = 0.
431            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  
432    
433            DO j = 1, knon            DO j = 1, knon
434               i = ni(j)               i = ni(j)
# Line 656  contains Line 437  contains
437                  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)
438                  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)
439                  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)  
440               END DO               END DO
441            END DO            END DO
442    
443            !cc diagnostic t, q a 2m et u, v a 10m            forall (k = 2:klev) coefh(ni(:knon), k) &
444                   = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypctsrf(:knon)
445    
446              ! diagnostic t, q a 2m et u, v a 10m
447    
448            DO j = 1, knon            DO j = 1, knon
449               i = ni(j)               i = ni(j)
450               uzon(j) = yu(j, 1) + y_d_u(j, 1)               u1(j) = yu(j, 1) + y_d_u(j, 1)
451               vmer(j) = yv(j, 1) + y_d_v(j, 1)               v1(j) = yv(j, 1) + y_d_v(j, 1)
452               tair1(j) = yt(j, 1) + y_d_t(j, 1)               tair1(j) = yt(j, 1) + y_d_t(j, 1)
453               qair1(j) = yq(j, 1) + y_d_q(j, 1)               qair1(j) = yq(j, 1) + y_d_q(j, 1)
454               zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &               zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
455                    1)))*(ypaprs(j, 1)-ypplay(j, 1))                    1))) * (ypaprs(j, 1)-ypplay(j, 1))
456               tairsol(j) = yts(j) + y_d_ts(j)               tairsol(j) = yts(j) + y_d_ts(j)
457               rugo1(j) = yrugos(j)               rugo1(j) = yrugos(j)
458               IF (nsrf == is_oce) THEN               IF (nsrf == is_oce) THEN
459                  rugo1(j) = rugos(i, nsrf)                  rugo1(j) = frugs(i, nsrf)
460               END IF               END IF
461               psfce(j) = ypaprs(j, 1)               psfce(j) = ypaprs(j, 1)
462               patm(j) = ypplay(j, 1)               patm(j) = ypplay(j, 1)
   
              qairsol(j) = yqsurf(j)  
463            END DO            END DO
464    
465            CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &            CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), qair1, &
466                 zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &                 zgeo1, tairsol, yqsurf(:knon), rugo1, psfce, patm, yt2m, yq2m, &
467                 yt10m, yq10m, yu10m, yustar)                 yt10m, yq10m, wind10m(:knon), ustar(:knon))
468    
469            DO j = 1, knon            DO j = 1, knon
470               i = ni(j)               i = ni(j)
471               t2m(i, nsrf) = yt2m(j)               t2m(i, nsrf) = yt2m(j)
472               q2m(i, nsrf) = yq2m(j)               q2m(i, nsrf) = yq2m(j)
473    
474               ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman               u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
475               u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)                    / sqrt(u1(j)**2 + v1(j)**2)
476               v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)               v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
477                      / sqrt(u1(j)**2 + v1(j)**2)
478            END DO            END DO
479    
480            CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
481                 y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &                 y_flux_q(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
482                 ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)                 yq(:knon, :), ypblh(:knon), ycapcl, yoliqcl, ycteicl, ypblt, &
483                   ytherm, ylcl)
484    
485            DO j = 1, knon            DO j = 1, knon
486               i = ni(j)               i = ni(j)
# Line 709  contains Line 491  contains
491               cteicl(i, nsrf) = ycteicl(j)               cteicl(i, nsrf) = ycteicl(j)
492               pblt(i, nsrf) = ypblt(j)               pblt(i, nsrf) = ypblt(j)
493               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)  
494            END DO            END DO
495    
496            DO j = 1, knon            IF (iflag_pbl >= 6) q2(ni(:knon), :, nsrf) = yq2(:knon, :)
497               DO k = 1, klev + 1         else
498                  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  
499         end IF if_knon         end IF if_knon
500      END DO loop_surface      END DO loop_surface
501    
502      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
503        frugs(:, is_oce) = rugmer
504        pctsrf(:, is_oce) = pctsrf_new_oce
505        pctsrf(:, is_sic) = pctsrf_new_sic
506    
507      rugos(:, is_oce) = rugmer      CALL histwrite_phy("run_off_lic", run_off_lic)
     pctsrf = pctsrf_new  
508    
509    END SUBROUTINE clmain    END SUBROUTINE pbl_surface
510    
511  end module clmain_m  end module pbl_surface_m
Legend:
Removed from v.72  
changed lines
  Added in v.310
  ViewVC Help
Powered by ViewVC 1.1.21