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