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

Legend:
Removed from v.61  
changed lines
  Added in v.328

  ViewVC Help
Powered by ViewVC 1.1.21