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

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

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

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