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

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

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

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