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

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

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

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