/[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 61 by guez, Fri Apr 20 14:58:43 2012 UTC trunk/phylmd/Interface_surf/pbl_surface.f revision 303 by guez, Thu Sep 6 14:25:07 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, pplay, fsnow, qsurf, evap, falbe, fluxlat, &
9         soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil,&         rain_fall, snow_f, fsolsw, fsollw, frugs, agesno, rugoro, d_t, d_q, &
10         qsol, paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat,&         d_u, d_v, d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, &
11         rain_f, snow_f, solsw, sollw, sollwdown, fder, rlon, rlat, cufi,&         dflux_t, dflux_q, coefh, t2m, q2m, u10m_srf, v10m_srf, pblh, capcl, &
12         cvfi, rugos, debut, lafin, agesno, rugoro, d_t, d_q, d_u, d_v,&         oliqcl, cteicl, pblt, therm, plcl, fqcalving, ffonte, run_off_lic_0)
13         d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2,&  
14         dflux_t, dflux_q, zcoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh,&      ! From phylmd/clmain.F, version 1.6, 2005/11/16 14:47:19
15         capcl, oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl,&      ! Author: Z. X. Li (LMD/CNRS)
16         fqcalving, ffonte, run_off_lic_0, flux_o, flux_g, tslab, seaice)      ! Date: Aug. 18th, 1993
17        ! Objet : interface de couche limite (diffusion verticale)
18      ! From phylmd/clmain.F, version 1.6 2005/11/16 14:47:19  
19      ! Author: Z.X. Li (LMD/CNRS), date: 1993/08/18      ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
20      ! Objet : interface de "couche limite" (diffusion verticale)      ! de la couche limite pour les traceurs se fait avec "cltrac" et
21        ! ne tient pas compte de la diff\'erentiation des sous-fractions
22      ! Tout ce qui a trait aux traceurs est dans "phytrac" maintenant.      ! de sol.
     ! Pour l'instant le calcul de la couche limite pour les traceurs  
     ! se fait avec "cltrac" et ne tient pas compte de la différentiation  
     ! des sous-fractions de sol.  
   
     ! Pour pouvoir extraire les coefficients d'échanges et le vent  
     ! dans la première couche, trois champs supplémentaires ont été  
     ! créés : "zcoefh", "zu1" et "zv1". Pour l'instant nous avons  
     ! moyenné les valeurs de ces trois champs sur les 4 sous-surfaces  
     ! du modèle. Dans l'avenir, si les informations des sous-surfaces  
     ! doivent être prises en compte, il faudra sortir ces mêmes champs  
     ! en leur ajoutant une dimension, c'est-à-dire "nbsrf" (nombre de  
     ! sous-surfaces).  
23    
24      use calendar, ONLY : ymds2ju      use cdrag_m, only: cdrag
25      use clqh_m, only: clqh      use clqh_m, only: clqh
26      use coefkz_m, only: coefkz      use clvent_m, only: clvent
27      use coefkzmin_m, only: coefkzmin      use coef_diff_turb_m, only: coef_diff_turb
28      USE conf_phys_m, ONLY : iflag_pbl      USE conf_gcm_m, ONLY: lmt_pas
29      USE dimens_m, ONLY : iim, jjm      USE conf_phys_m, ONLY: iflag_pbl
30      USE dimphy, ONLY : klev, klon, zmasq      USE dimphy, ONLY: klev, klon
31      USE dimsoil, ONLY : nsoilmx      USE dimsoil, ONLY: nsoilmx
     USE dynetat0_m, ONLY : day_ini  
     USE gath_cpl, ONLY : gath2cpl  
32      use hbtm_m, only: hbtm      use hbtm_m, only: hbtm
33      USE histsync_m, ONLY : histsync      USE histwrite_phy_m, ONLY: histwrite_phy
34      USE histbeg_totreg_m, ONLY : histbeg_totreg      USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
35      USE histend_m, ONLY : histend      USE interfoce_lim_m, ONLY: interfoce_lim
36      USE histdef_m, ONLY : histdef      use phyetat0_m, only: zmasq
37      use histwrite_m, only: histwrite      use stdlevvar_m, only: stdlevvar
38      USE indicesol, ONLY : epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf      USE suphec_m, ONLY: rd, rg
39      USE conf_gcm_m, ONLY : prt_level      use time_phylmdz, only: itap
40      USE suphec_m, ONLY : rd, rg, rkappa  
41      USE temps, ONLY : annee_ref, itau_phy      REAL, INTENT(inout):: pctsrf(klon, nbsrf)
42      use yamada4_m, only: yamada4      ! tableau des pourcentages de surface de chaque maille
43    
44      ! Arguments:      REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
45        REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg / kg)
46      REAL, INTENT (IN) :: dtime ! interval du temps (secondes)      REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
47      REAL date0      INTEGER, INTENT(IN):: julien ! jour de l'annee en cours
48      ! date0----input-R- jour initial      REAL, intent(in):: mu0(klon) ! cosinus de l'angle solaire zenithal    
49      INTEGER, INTENT (IN) :: itap      REAL, INTENT(IN):: ftsol(:, :) ! (klon, nbsrf) temp\'erature du sol (en K)
50      ! itap-----input-I- numero du pas de temps      REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
51      REAL, INTENT(IN):: t(klon, klev), q(klon, klev)  
52      ! t--------input-R- temperature (K)      REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
53      ! q--------input-R- vapeur d'eau (kg/kg)      ! soil temperature of surface fraction
54      REAL, INTENT (IN):: u(klon, klev), v(klon, klev)  
55      ! u--------input-R- vitesse u      REAL, INTENT(inout):: qsol(:) ! (klon)
56      ! v--------input-R- vitesse v      ! column-density of water in soil, in kg m-2
57      REAL, INTENT (IN):: paprs(klon, klev+1)  
58      ! paprs----input-R- pression a intercouche (Pa)      REAL, INTENT(IN):: paprs(klon, klev + 1) ! pression a intercouche (Pa)
59      REAL, INTENT (IN):: pplay(klon, klev)      REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
60      ! pplay----input-R- pression au milieu de couche (Pa)      REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse
     REAL, INTENT (IN):: rlon(klon), rlat(klon)  
     ! rlat-----input-R- latitude en degree  
     REAL cufi(klon), cvfi(klon)  
     ! cufi-----input-R- resolution des mailles en x (m)  
     ! cvfi-----input-R- resolution des mailles en y (m)  
     REAL d_t(klon, klev), d_q(klon, klev)  
     ! d_t------output-R- le changement pour "t"  
     ! d_q------output-R- le changement pour "q"  
     REAL d_u(klon, klev), d_v(klon, klev)  
     ! d_u------output-R- le changement pour "u"  
     ! d_v------output-R- le changement pour "v"  
     REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)  
     ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)  
     !                    (orientation positive vers le bas)  
     ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)  
     REAL dflux_t(klon), dflux_q(klon)  
     ! dflux_t derive du flux sensible  
     ! dflux_q derive du flux latent  
     !IM "slab" ocean  
     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 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)  
     REAL cdragh(klon), cdragm(klon)  
     ! jour de l'annee en cours                  
     INTEGER jour  
     REAL rmu0(klon) ! cosinus de l'angle solaire zenithal      
     ! taux CO2 atmosphere                      
     REAL co2_ppm  
     LOGICAL, INTENT (IN) :: debut  
     LOGICAL, INTENT (IN) :: lafin  
     LOGICAL ok_veget  
     CHARACTER (len=*), INTENT (IN) :: ocean  
     INTEGER npas, nexca  
   
     REAL pctsrf(klon, nbsrf)  
     REAL ts(klon, nbsrf)  
     ! ts-------input-R- temperature du sol (en Kelvin)  
     REAL d_ts(klon, nbsrf)  
     ! d_ts-----output-R- le changement pour "ts"  
     REAL snow(klon, nbsrf)  
61      REAL qsurf(klon, nbsrf)      REAL qsurf(klon, nbsrf)
62      REAL evap(klon, nbsrf)      REAL evap(klon, nbsrf)
63      REAL albe(klon, nbsrf)      REAL, intent(inout):: falbe(klon, nbsrf)
64      REAL alblw(klon, nbsrf)      REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)
65    
66      REAL fluxlat(klon, nbsrf)      REAL, intent(in):: rain_fall(klon)
67        ! liquid water mass flux (kg / m2 / s), positive down
68    
69      REAL rain_f(klon), snow_f(klon)      REAL, intent(in):: snow_f(klon)
70      REAL fder(klon)      ! solid water mass flux (kg / m2 / s), positive down
71    
72      REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)      REAL, INTENT(IN):: fsolsw(klon, nbsrf), fsollw(klon, nbsrf)
73      REAL rugos(klon, nbsrf)      REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)
74      ! rugos----input-R- longeur de rugosite (en m)      real agesno(klon, nbsrf)
75      ! la nouvelle repartition des surfaces sortie de l'interface      REAL, INTENT(IN):: rugoro(klon)
     REAL pctsrf_new(klon, nbsrf)  
76    
77      REAL zcoefh(klon, klev)      REAL, intent(out):: d_t(:, :), d_q(:, :) ! (klon, klev)
78      REAL zu1(klon)      ! changement pour t et q
     REAL zv1(klon)  
79    
80      !$$$ PB ajout pour soil      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
81      LOGICAL, INTENT (IN) :: soil_model      ! changement pour "u" et "v"
     !IM ajout seuils cdrm, cdrh  
     REAL cdmmax, cdhmax  
82    
83      REAL ksta, ksta_ter      REAL, intent(out):: d_ts(:, :) ! (klon, nbsrf) variation of ftsol
     LOGICAL ok_kzmin  
84    
85      REAL ftsoil(klon, nsoilmx, nbsrf)      REAL, intent(out):: flux_t(klon, nbsrf)
86      REAL ytsoil(klon, nsoilmx)      ! flux de chaleur sensible (c_p T) (W / m2) (orientation positive
87      REAL qsol(klon)      ! vers le bas) à la surface
88    
89        REAL, intent(out):: flux_q(klon, nbsrf)
90        ! flux de vapeur d'eau (kg / m2 / s) à la surface
91    
92        REAL, intent(out):: flux_u(klon, nbsrf), flux_v(klon, nbsrf)
93        ! tension du vent (flux turbulent de vent) à la surface, en Pa
94    
95        REAL, INTENT(out):: cdragh(klon), cdragm(klon)
96        real q2(klon, klev + 1, nbsrf)
97    
98        ! Ocean slab:
99        REAL, INTENT(out):: dflux_t(klon) ! derive du flux sensible
100        REAL, INTENT(out):: dflux_q(klon) ! derive du flux latent
101    
102        REAL, intent(out):: coefh(:, 2:) ! (klon, 2:klev)
103        ! Pour pouvoir extraire les coefficients d'\'echange, le champ
104        ! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
105        ! ce champ sur les quatre sous-surfaces du mod\`ele.
106    
107        REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)
108    
109        REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)
110        ! composantes du vent \`a 10m sans spirale d'Ekman
111    
112        ! Ionela Musat. Cf. Anne Mathieu : planetary boundary layer, hbtm.
113        ! Comme les autres diagnostics on cumule dans physiq ce qui permet
114        ! de sortir les grandeurs par sous-surface.
115        REAL pblh(klon, nbsrf) ! height of planetary boundary layer
116        REAL capcl(klon, nbsrf)
117        REAL oliqcl(klon, nbsrf)
118        REAL cteicl(klon, nbsrf)
119        REAL, INTENT(inout):: pblt(klon, nbsrf) ! T au nveau HCL
120        REAL therm(klon, nbsrf)
121        REAL plcl(klon, nbsrf)
122    
123        REAL, intent(out):: fqcalving(klon, nbsrf)
124        ! flux d'eau "perdue" par la surface et necessaire pour limiter la
125        ! hauteur de neige, en kg / m2 / s
126    
127      EXTERNAL clvent, calbeta, cltrac      real ffonte(klon, nbsrf) ! flux thermique utilise pour fondre la neige
128        REAL, intent(inout):: run_off_lic_0(:) ! (klon)
129    
130      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      ! Local:
131    
132        ! la nouvelle repartition des surfaces sortie de l'interface
133        REAL, save:: pctsrf_new_oce(klon)
134        REAL, save:: pctsrf_new_sic(klon)
135    
136        REAL y_fqcalving(klon), y_ffonte(klon)
137        real y_run_off_lic_0(klon), y_run_off_lic(klon)
138        REAL run_off_lic(klon) ! ruissellement total
139        REAL rugmer(klon)
140        REAL ytsoil(klon, nsoilmx)
141        REAL yts(klon), ypct(klon), yz0_new(klon)
142        real yrugos(klon) ! longueur de rugosite (en m)
143      REAL yalb(klon)      REAL yalb(klon)
144      REAL yalblw(klon)      REAL snow(klon), yqsurf(klon), yagesno(klon)
145      REAL yu1(klon), yv1(klon)      real yqsol(klon) ! column-density of water in soil, in kg m-2
146      ! on rajoute en output yu1 et yv1 qui sont les vents dans      REAL yrain_f(klon) ! liquid water mass flux (kg / m2 / s), positive down
147      ! la premiere couche      REAL ysnow_f(klon) ! solid water mass flux (kg / m2 / s), positive down
     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)  
148      REAL yrugm(klon), yrads(klon), yrugoro(klon)      REAL yrugm(klon), yrads(klon), yrugoro(klon)
   
149      REAL yfluxlat(klon)      REAL yfluxlat(klon)
   
150      REAL y_d_ts(klon)      REAL y_d_ts(klon)
151      REAL y_d_t(klon, klev), y_d_q(klon, klev)      REAL y_d_t(klon, klev), y_d_q(klon, klev)
152      REAL y_d_u(klon, klev), y_d_v(klon, klev)      REAL y_d_u(klon, klev), y_d_v(klon, klev)
153      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)      REAL y_flux_t(klon), y_flux_q(klon)
154      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)      REAL y_flux_u(klon), y_flux_v(klon)
155      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
156      REAL ycoefh(klon, klev), ycoefm(klon, klev)      REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
157        real ycdragh(klon), ycdragm(klon)
158      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
159      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
160      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)      REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
161        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)  
162      REAL delp(klon, klev)      REAL delp(klon, klev)
163      INTEGER i, k, nsrf      INTEGER i, k, nsrf
   
164      INTEGER ni(klon), knon, j      INTEGER ni(klon), knon, j
165    
166      REAL pctsrf_pot(klon, nbsrf)      REAL pctsrf_pot(klon, nbsrf)
167      ! "pourcentage potentiel" pour tenir compte des éventuelles      ! "pourcentage potentiel" pour tenir compte des \'eventuelles
168      ! apparitions ou disparitions de la glace de mer      ! apparitions ou disparitions de la glace de mer
169    
170      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.      REAL yt2m(klon), yq2m(klon), wind10m(klon)
171        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  
172    
173      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  
174      REAL ypblh(klon)      REAL ypblh(klon)
175      REAL ylcl(klon)      REAL ylcl(klon)
176      REAL ycapcl(klon)      REAL ycapcl(klon)
# Line 265  contains Line 178  contains
178      REAL ycteicl(klon)      REAL ycteicl(klon)
179      REAL ypblt(klon)      REAL ypblt(klon)
180      REAL ytherm(klon)      REAL ytherm(klon)
181      REAL ytrmb1(klon)      REAL u1(klon), v1(klon)
     REAL ytrmb2(klon)  
     REAL ytrmb3(klon)  
     REAL y_cd_h(klon), y_cd_m(klon)  
     REAL uzon(klon), vmer(klon)  
182      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
183      REAL psfce(klon), patm(klon)      REAL psfce(klon), patm(klon)
184    
185      REAL qairsol(klon), zgeo1(klon)      REAL qairsol(klon), zgeo1(klon)
186      REAL rugo1(klon)      REAL rugo1(klon)
187        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'  
188    
189      !------------------------------------------------------------      !------------------------------------------------------------
190    
191      ytherm = 0.      ytherm = 0.
192    
     IF (debugindex .AND. first_appel) THEN  
        first_appel = .FALSE.  
   
        ! initialisation sorties netcdf  
   
        idayref = day_ini  
        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  
   
193      DO k = 1, klev ! epaisseur de couche      DO k = 1, klev ! epaisseur de couche
194         DO i = 1, klon         DO i = 1, klon
195            delp(i, k) = paprs(i, k) - paprs(i, k+1)            delp(i, k) = paprs(i, k) - paprs(i, k + 1)
196         END DO         END DO
197      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  
198    
199      ! Initialization:      ! Initialization:
200      rugmer = 0.      rugmer = 0.
# Line 337  contains Line 202  contains
202      cdragm = 0.      cdragm = 0.
203      dflux_t = 0.      dflux_t = 0.
204      dflux_q = 0.      dflux_q = 0.
     zu1 = 0.  
     zv1 = 0.  
205      ypct = 0.      ypct = 0.
     yts = 0.  
     ysnow = 0.  
206      yqsurf = 0.      yqsurf = 0.
     yalb = 0.  
     yalblw = 0.  
207      yrain_f = 0.      yrain_f = 0.
208      ysnow_f = 0.      ysnow_f = 0.
     yfder = 0.  
     ytaux = 0.  
     ytauy = 0.  
     ysolsw = 0.  
     ysollw = 0.  
     ysollwdown = 0.  
209      yrugos = 0.      yrugos = 0.
     yu1 = 0.  
     yv1 = 0.  
     yrads = 0.  
210      ypaprs = 0.      ypaprs = 0.
211      ypplay = 0.      ypplay = 0.
212      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.  
213      yrugoro = 0.      yrugoro = 0.
     ! -- LOOP  
     yu10mx = 0.  
     yu10my = 0.  
     ywindsp = 0.  
     ! -- LOOP  
214      d_ts = 0.      d_ts = 0.
     !§§§ PB  
     yfluxlat = 0.  
215      flux_t = 0.      flux_t = 0.
216      flux_q = 0.      flux_q = 0.
217      flux_u = 0.      flux_u = 0.
218      flux_v = 0.      flux_v = 0.
219        fluxlat = 0.
220      d_t = 0.      d_t = 0.
221      d_q = 0.      d_q = 0.
222      d_u = 0.      d_u = 0.
223      d_v = 0.      d_v = 0.
224      zcoefh = 0.      coefh = 0.
225        fqcalving = 0.
226      ! Boucler sur toutes les sous-fractions du sol:      run_off_lic = 0.
227    
228      ! Initialisation des "pourcentages potentiels". On considère ici qu'on      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
229      ! peut avoir potentiellement de la glace sur tout le domaine océanique      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
230      ! (à affiner)      ! (\`a affiner).
231    
232      pctsrf_pot = pctsrf      pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
233        pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
234      pctsrf_pot(:, is_oce) = 1. - zmasq      pctsrf_pot(:, is_oce) = 1. - zmasq
235      pctsrf_pot(:, is_sic) = 1. - zmasq      pctsrf_pot(:, is_sic) = 1. - zmasq
236    
237        ! Tester si c'est le moment de lire le fichier:
238        if (mod(itap - 1, lmt_pas) == 0) then
239           CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
240        endif
241    
242        ! Boucler sur toutes les sous-fractions du sol:
243    
244      loop_surface: DO nsrf = 1, nbsrf      loop_surface: DO nsrf = 1, nbsrf
245         ! Chercher les indices :         ! Define ni and knon:
246          
247         ni = 0         ni = 0
248         knon = 0         knon = 0
249    
250         DO i = 1, klon         DO i = 1, klon
251            ! Pour déterminer le domaine à traiter, on utilise les surfaces            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
252            ! "potentielles"            ! "potentielles"
253            IF (pctsrf_pot(i, nsrf) > epsfra) THEN            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
254               knon = knon + 1               knon = knon + 1
# Line 413  contains Line 256  contains
256            END IF            END IF
257         END DO         END DO
258    
259         ! variables pour avoir une sortie IOIPSL des INDEX         if_knon: IF (knon /= 0) then
        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  
   
        IF (knon == 0) CYCLE  
   
        DO j = 1, knon  
           i = ni(j)  
           ypct(j) = pctsrf(i, nsrf)  
           yts(j) = ts(i, nsrf)  
           ytslab(i) = tslab(i)  
           ysnow(j) = snow(i, nsrf)  
           yqsurf(j) = qsurf(i, nsrf)  
           yalb(j) = albe(i, nsrf)  
           yalblw(j) = alblw(i, nsrf)  
           yrain_f(j) = rain_f(i)  
           ysnow_f(j) = snow_f(i)  
           yagesno(j) = agesno(i, nsrf)  
           yfder(j) = fder(i)  
           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  
   
        ! IF bucket model for continent, copy soil water content  
        IF (nsrf == is_ter .AND. .NOT. ok_veget) THEN  
           DO j = 1, knon  
              i = ni(j)  
              yqsol(j) = qsol(i)  
           END DO  
        ELSE  
           yqsol = 0.  
        END IF  
        !$$$ PB ajour pour soil  
        DO k = 1, nsoilmx  
           DO j = 1, knon  
              i = ni(j)  
              ytsoil(j, k) = ftsoil(i, k, nsrf)  
           END DO  
        END DO  
        DO k = 1, klev  
260            DO j = 1, knon            DO j = 1, knon
261               i = ni(j)               i = ni(j)
262               ypaprs(j, k) = paprs(i, k)               ypct(j) = pctsrf(i, nsrf)
263               ypplay(j, k) = pplay(i, k)               yts(j) = ftsol(i, nsrf)
264               ydelp(j, k) = delp(i, k)               snow(j) = fsnow(i, nsrf)
265               yu(j, k) = u(i, k)               yqsurf(j) = qsurf(i, nsrf)
266               yv(j, k) = v(i, k)               yalb(j) = falbe(i, nsrf)
267               yt(j, k) = t(i, k)               yrain_f(j) = rain_fall(i)
268               yq(j, k) = q(i, k)               ysnow_f(j) = snow_f(i)
269            END DO               yagesno(j) = agesno(i, nsrf)
270         END DO               yrugos(j) = frugs(i, nsrf)
271                 yrugoro(j) = rugoro(i)
272         ! calculer Cdrag et les coefficients d'echange               yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)
273         CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,&               ypaprs(j, klev + 1) = paprs(i, klev + 1)
274              yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)               y_run_off_lic_0(j) = run_off_lic_0(i)
        IF (iflag_pbl == 1) THEN  
           CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)  
           DO k = 1, klev  
              DO i = 1, knon  
                 ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))  
                 ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))  
              END DO  
           END DO  
        END IF  
   
        ! on seuille ycoefm et ycoefh  
        IF (nsrf == is_oce) THEN  
           DO j = 1, knon  
              ycoefm(j, 1) = min(ycoefm(j, 1), cdmmax)  
              ycoefh(j, 1) = min(ycoefh(j, 1), cdhmax)  
275            END DO            END DO
        END IF  
276    
277         IF (ok_kzmin) THEN            ! For continent, copy soil water content
278            ! Calcul d'une diffusion minimale pour les conditions tres stables            IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
           CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, ycoefm(:, 1), &  
                ycoefm0, ycoefh0)  
279    
280            DO k = 1, klev            ytsoil(:knon, :) = ftsoil(ni(:knon), :, nsrf)
              DO i = 1, knon  
                 ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))  
                 ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))  
              END DO  
           END DO  
        END IF  
281    
        IF (iflag_pbl >= 3) THEN  
           ! MELLOR ET YAMADA adapté à Mars, Richard Fournier et Frédéric Hourdin  
           yzlay(1:knon, 1) = rd*yt(1:knon, 1)/(0.5*(ypaprs(1:knon, &  
                1)+ypplay(1:knon, 1)))*(ypaprs(1:knon, 1)-ypplay(1: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  
282            DO k = 1, klev            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(1:knon, klev+1) = 2.*yzlay(1:knon, klev) - yzlay(1: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  
283               DO j = 1, knon               DO j = 1, knon
284                  i = ni(j)                  i = ni(j)
285                  yq2(j, k) = q2(i, k, nsrf)                  ypaprs(j, k) = paprs(i, k)
286                    ypplay(j, k) = pplay(i, k)
287                    ydelp(j, k) = delp(i, k)
288                    yu(j, k) = u(i, k)
289                    yv(j, k) = v(i, k)
290                    yt(j, k) = t(i, k)
291                    yq(j, k) = q(i, k)
292               END DO               END DO
293            END DO            END DO
294    
295            y_cd_m(1:knon) = ycoefm(1:knon, 1)            ! Calculer les géopotentiels de chaque couche:
           y_cd_h(1:knon) = ycoefh(1:knon, 1)  
           CALL ustarhb(knon, yu, yv, y_cd_m, yustar)  
296    
297            IF (prt_level>9) THEN            zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
298               PRINT *, 'USTAR = ', yustar                 + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
           END IF  
299    
300            ! iflag_pbl peut être utilisé comme longueur de mélange            DO k = 2, klev
301                 zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
302                      * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
303                      * (ypplay(:knon, k - 1) - ypplay(:knon, k))
304              ENDDO
305    
306              CALL cdrag(nsrf, sqrt(yu(:knon, 1)**2 + yv(:knon, 1)**2), &
307                   yt(:knon, 1), yq(:knon, 1), zgeop(:knon, 1), ypaprs(:knon, 1), &
308                   yts(:knon), yqsurf(:knon), yrugos(:knon), ycdragm(:knon), &
309                   ycdragh(:knon))
310    
311              IF (iflag_pbl == 1) THEN
312                 ycdragm(:knon) = max(ycdragm(:knon), 0.)
313                 ycdragh(:knon) = max(ycdragh(:knon), 0.)
314              end IF
315    
316            IF (iflag_pbl >= 11) THEN            ! on met un seuil pour ycdragm et ycdragh
317               CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &            IF (nsrf == is_oce) THEN
318                    yu, yv, yteta, y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, &               ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
319                    iflag_pbl)               ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
           ELSE  
              CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &  
                   y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)  
320            END IF            END IF
321    
322            ycoefm(1:knon, 1) = y_cd_m(1:knon)            IF (iflag_pbl >= 6) yq2(:knon, :) = q2(ni(:knon), :, nsrf)
323            ycoefh(1:knon, 1) = y_cd_h(1:knon)            call coef_diff_turb(nsrf, ni(:knon), ypaprs(:knon, :), &
324            ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)                 ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
325            ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)                 yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
326         END IF                 ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
327              
328         ! calculer la diffusion des vitesses "u" et "v"            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
329         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &                 ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
330              ydelp, y_d_u, y_flux_u)                 ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
331         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &                 y_flux_u(:knon))
332              ydelp, y_d_v, y_flux_v)            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
333                   ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
334         ! pour le couplage                 ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
335         ytaux = y_flux_u(:, 1)                 y_flux_v(:knon))
336         ytauy = y_flux_v(:, 1)  
337              CALL clqh(julien, nsrf, ni(:knon), ytsoil(:knon, :), yqsol(:knon), &
338         ! calculer la diffusion de "q" et de "h"                 mu0(ni(:knon)), yrugos(:knon), yrugoro(:knon), yu(:knon, 1), &
339         CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,&                 yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), yt(:knon, :), &
340              cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil,&                 yq(:knon, :), yts(:knon), ypaprs(:knon, :), ypplay(:knon, :), &
341              yqsol, ok_veget, ocean, npas, nexca, rmu0, co2_ppm, yrugos,&                 ydelp(:knon, :), yrads(:knon), yalb(:knon), snow(:knon), &
342              yrugoro, yu1, yv1, ycoefh, yt, yq, yts, ypaprs, ypplay,&                 yqsurf(:knon), yrain_f(:knon), ysnow_f(:knon), yfluxlat(:knon), &
343              ydelp, yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, &                 pctsrf_new_sic(ni(:knon)), yagesno(:knon), y_d_t(:knon, :), &
344              yfder, ytaux, ytauy, ywindsp, ysollw, ysollwdown, ysolsw,&                 y_d_q(:knon, :), y_d_ts(:knon), yz0_new(:knon), &
345              yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts,&                 y_flux_t(:knon), y_flux_q(:knon), y_dflux_t(:knon), &
346              yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,&                 y_dflux_q(:knon), y_fqcalving(:knon), y_ffonte(:knon), &
347              y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g,&                 y_run_off_lic_0(:knon), y_run_off_lic(:knon))
             ytslab, y_seaice)  
   
        ! calculer la longueur de rugosite sur ocean  
        yrugm = 0.  
        IF (nsrf == is_oce) THEN  
           DO j = 1, knon  
              yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &  
                   0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2))  
              yrugm(j) = max(1.5E-05, yrugm(j))  
           END DO  
        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  
348    
349         DO k = 1, klev            ! calculer la longueur de rugosite sur ocean
           DO j = 1, knon  
              i = ni(j)  
              ycoefh(j, k) = ycoefh(j, k)*ypct(j)  
              ycoefm(j, k) = ycoefm(j, k)*ypct(j)  
              y_d_t(j, k) = y_d_t(j, k)*ypct(j)  
              y_d_q(j, k) = y_d_q(j, k)*ypct(j)  
              flux_t(i, k, nsrf) = y_flux_t(j, k)  
              flux_q(i, k, nsrf) = y_flux_q(j, k)  
              flux_u(i, k, nsrf) = y_flux_u(j, k)  
              flux_v(i, k, nsrf) = y_flux_v(j, k)  
              y_d_u(j, k) = y_d_u(j, k)*ypct(j)  
              y_d_v(j, k) = y_d_v(j, k)*ypct(j)  
           END DO  
        END DO  
350    
351         evap(:, nsrf) = -flux_q(:, 1, nsrf)            yrugm = 0.
352    
        albe(:, nsrf) = 0.  
        alblw(:, nsrf) = 0.  
        snow(:, nsrf) = 0.  
        qsurf(:, nsrf) = 0.  
        rugos(:, nsrf) = 0.  
        fluxlat(:, nsrf) = 0.  
        DO j = 1, knon  
           i = ni(j)  
           d_ts(i, nsrf) = y_d_ts(j)  
           albe(i, nsrf) = yalb(j)  
           alblw(i, nsrf) = yalblw(j)  
           snow(i, nsrf) = ysnow(j)  
           qsurf(i, nsrf) = yqsurf(j)  
           rugos(i, nsrf) = yz0_new(j)  
           fluxlat(i, nsrf) = yfluxlat(j)  
353            IF (nsrf == is_oce) THEN            IF (nsrf == is_oce) THEN
354               rugmer(i) = yrugm(j)               DO j = 1, knon
355               rugos(i, nsrf) = yrugm(j)                  yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
356                         / rg + 0.11 * 14E-6 &
357                         / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
358                    yrugm(j) = max(1.5E-05, yrugm(j))
359                 END DO
360            END IF            END IF
           agesno(i, nsrf) = yagesno(j)  
           fqcalving(i, nsrf) = y_fqcalving(j)  
           ffonte(i, nsrf) = y_ffonte(j)  
           cdragh(i) = cdragh(i) + ycoefh(j, 1)  
           cdragm(i) = cdragm(i) + ycoefm(j, 1)  
           dflux_t(i) = dflux_t(i) + y_dflux_t(j)  
           dflux_q(i) = dflux_q(i) + y_dflux_q(j)  
           zu1(i) = zu1(i) + yu1(j)  
           zv1(i) = zv1(i) + yv1(j)  
        END DO  
        IF (nsrf == is_ter) THEN  
           DO j = 1, knon  
              i = ni(j)  
              qsol(i) = yqsol(j)  
           END DO  
        END IF  
        IF (nsrf == is_lic) THEN  
           DO j = 1, knon  
              i = ni(j)  
              run_off_lic_0(i) = y_run_off_lic_0(j)  
           END DO  
        END IF  
        !$$$ PB ajout pour soil  
        ftsoil(:, :, nsrf) = 0.  
        DO k = 1, nsoilmx  
           DO j = 1, knon  
              i = ni(j)  
              ftsoil(i, k, nsrf) = ytsoil(j, k)  
           END DO  
        END DO  
361    
        DO j = 1, knon  
           i = ni(j)  
362            DO k = 1, klev            DO k = 1, klev
363               d_t(i, k) = d_t(i, k) + y_d_t(j, k)               DO j = 1, knon
364               d_q(i, k) = d_q(i, k) + y_d_q(j, k)                  i = ni(j)
365               d_u(i, k) = d_u(i, k) + y_d_u(j, k)                  y_d_t(j, k) = y_d_t(j, k) * ypct(j)
366               d_v(i, k) = d_v(i, k) + y_d_v(j, k)                  y_d_q(j, k) = y_d_q(j, k) * ypct(j)
367               zcoefh(i, k) = zcoefh(i, k) + ycoefh(j, k)                  y_d_u(j, k) = y_d_u(j, k) * ypct(j)
368                    y_d_v(j, k) = y_d_v(j, k) * ypct(j)
369                 END DO
370            END DO            END DO
        END DO  
371    
372         !cc diagnostic t, q a 2m et u, v a 10m            flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
373              flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
374         DO j = 1, knon            flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
375            i = ni(j)            flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
376            uzon(j) = yu(j, 1) + y_d_u(j, 1)  
377            vmer(j) = yv(j, 1) + y_d_v(j, 1)            evap(:, nsrf) = -flux_q(:, nsrf)
378            tair1(j) = yt(j, 1) + y_d_t(j, 1)  
379            qair1(j) = yq(j, 1) + y_d_q(j, 1)            falbe(:, nsrf) = 0.
380            zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &            fsnow(:, nsrf) = 0.
381                 1)))*(ypaprs(j, 1)-ypplay(j, 1))            qsurf(:, nsrf) = 0.
382            tairsol(j) = yts(j) + y_d_ts(j)            frugs(:, nsrf) = 0.
383            rugo1(j) = yrugos(j)            DO j = 1, knon
384            IF (nsrf == is_oce) THEN               i = ni(j)
385               rugo1(j) = rugos(i, nsrf)               d_ts(i, nsrf) = y_d_ts(j)
386                 falbe(i, nsrf) = yalb(j)
387                 fsnow(i, nsrf) = snow(j)
388                 qsurf(i, nsrf) = yqsurf(j)
389                 frugs(i, nsrf) = yz0_new(j)
390                 fluxlat(i, nsrf) = yfluxlat(j)
391                 IF (nsrf == is_oce) THEN
392                    rugmer(i) = yrugm(j)
393                    frugs(i, nsrf) = yrugm(j)
394                 END IF
395                 agesno(i, nsrf) = yagesno(j)
396                 fqcalving(i, nsrf) = y_fqcalving(j)
397                 ffonte(i, nsrf) = y_ffonte(j)
398                 cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
399                 cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
400                 dflux_t(i) = dflux_t(i) + y_dflux_t(j) * ypct(j)
401                 dflux_q(i) = dflux_q(i) + y_dflux_q(j) * ypct(j)
402              END DO
403              IF (nsrf == is_ter) THEN
404                 qsol(ni(:knon)) = yqsol(:knon)
405              else IF (nsrf == is_lic) THEN
406                 DO j = 1, knon
407                    i = ni(j)
408                    run_off_lic_0(i) = y_run_off_lic_0(j)
409                    run_off_lic(i) = y_run_off_lic(j)
410                 END DO
411            END IF            END IF
           psfce(j) = ypaprs(j, 1)  
           patm(j) = ypplay(j, 1)  
412    
413            qairsol(j) = yqsurf(j)            ftsoil(:, :, nsrf) = 0.
414         END DO            ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)
415    
416         CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &            DO j = 1, knon
417              tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &               i = ni(j)
418              yu10m, yustar)               DO k = 1, klev
419                    d_t(i, k) = d_t(i, k) + y_d_t(j, k)
420         DO j = 1, knon                  d_q(i, k) = d_q(i, k) + y_d_q(j, k)
421            i = ni(j)                  d_u(i, k) = d_u(i, k) + y_d_u(j, k)
422            t2m(i, nsrf) = yt2m(j)                  d_v(i, k) = d_v(i, k) + y_d_v(j, k)
423            q2m(i, nsrf) = yq2m(j)               END DO
424              END DO
           ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman  
           u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)  
           v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)  
425    
426         END DO            forall (k = 2:klev) coefh(ni(:knon), k) &
427                   = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypct(:knon)
428    
429         DO i = 1, knon            ! diagnostic t, q a 2m et u, v a 10m
           y_cd_h(i) = ycoefh(i, 1)  
           y_cd_m(i) = ycoefm(i, 1)  
        END DO  
        CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &  
             y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &  
             ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)  
   
        DO j = 1, knon  
           i = ni(j)  
           pblh(i, nsrf) = ypblh(j)  
           plcl(i, nsrf) = ylcl(j)  
           capcl(i, nsrf) = ycapcl(j)  
           oliqcl(i, nsrf) = yoliqcl(j)  
           cteicl(i, nsrf) = ycteicl(j)  
           pblt(i, nsrf) = ypblt(j)  
           therm(i, nsrf) = ytherm(j)  
           trmb1(i, nsrf) = ytrmb1(j)  
           trmb2(i, nsrf) = ytrmb2(j)  
           trmb3(i, nsrf) = ytrmb3(j)  
        END DO  
430    
431         DO j = 1, knon            DO j = 1, knon
           DO k = 1, klev + 1  
432               i = ni(j)               i = ni(j)
433               q2(i, k, nsrf) = yq2(j, k)               u1(j) = yu(j, 1) + y_d_u(j, 1)
434                 v1(j) = yv(j, 1) + y_d_v(j, 1)
435                 tair1(j) = yt(j, 1) + y_d_t(j, 1)
436                 qair1(j) = yq(j, 1) + y_d_q(j, 1)
437                 zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
438                      1))) * (ypaprs(j, 1)-ypplay(j, 1))
439                 tairsol(j) = yts(j) + y_d_ts(j)
440                 rugo1(j) = yrugos(j)
441                 IF (nsrf == is_oce) THEN
442                    rugo1(j) = frugs(i, nsrf)
443                 END IF
444                 psfce(j) = ypaprs(j, 1)
445                 patm(j) = ypplay(j, 1)
446    
447                 qairsol(j) = yqsurf(j)
448            END DO            END DO
449         END DO  
450         !IM "slab" ocean            CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), qair1, &
451         IF (nsrf == is_oce) THEN                 zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, &
452                   yq10m, wind10m(:knon), ustar(:knon))
453    
454            DO j = 1, knon            DO j = 1, knon
              ! on projette sur la grille globale  
455               i = ni(j)               i = ni(j)
456               IF (pctsrf_new(i, is_oce)>epsfra) THEN               t2m(i, nsrf) = yt2m(j)
457                  flux_o(i) = y_flux_o(j)               q2m(i, nsrf) = yq2m(j)
458               ELSE  
459                  flux_o(i) = 0.               u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
460               END IF                    / sqrt(u1(j)**2 + v1(j)**2)
461                 v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
462                      / sqrt(u1(j)**2 + v1(j)**2)
463            END DO            END DO
        END IF  
464    
465         IF (nsrf == is_sic) THEN            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
466                   y_flux_q(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
467                   yq(:knon, :), ypblh(:knon), ycapcl, yoliqcl, ycteicl, ypblt, &
468                   ytherm, ylcl)
469    
470            DO j = 1, knon            DO j = 1, knon
471               i = ni(j)               i = ni(j)
472               ! On pondère lorsque l'on fait le bilan au sol :               pblh(i, nsrf) = ypblh(j)
473               IF (pctsrf_new(i, is_sic)>epsfra) THEN               plcl(i, nsrf) = ylcl(j)
474                  flux_g(i) = y_flux_g(j)               capcl(i, nsrf) = ycapcl(j)
475               ELSE               oliqcl(i, nsrf) = yoliqcl(j)
476                  flux_g(i) = 0.               cteicl(i, nsrf) = ycteicl(j)
477               END IF               pblt(i, nsrf) = ypblt(j)
478                 therm(i, nsrf) = ytherm(j)
479            END DO            END DO
480    
481         END IF            IF (iflag_pbl >= 6) q2(ni(:knon), :, nsrf) = yq2(:knon, :)
482         IF (ocean == 'slab  ') THEN         else
483            IF (nsrf == is_oce) THEN            fsnow(:, nsrf) = 0.
484               tslab(1:klon) = ytslab(1:klon)         end IF if_knon
              seaice(1:klon) = y_seaice(1:klon)  
           END IF  
        END IF  
485      END DO loop_surface      END DO loop_surface
486    
487      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
488        frugs(:, is_oce) = rugmer
489        pctsrf(:, is_oce) = pctsrf_new_oce
490        pctsrf(:, is_sic) = pctsrf_new_sic
491    
492      rugos(:, is_oce) = rugmer      CALL histwrite_phy("run_off_lic", run_off_lic)
     pctsrf = pctsrf_new  
493    
494    END SUBROUTINE clmain    END SUBROUTINE pbl_surface
495    
496  end module clmain_m  end module pbl_surface_m
Legend:
Removed from v.61  
changed lines
  Added in v.303
  ViewVC Help
Powered by ViewVC 1.1.21