/[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 38 by guez, Thu Jan 6 17:52:19 2011 UTC trunk/phylmd/Interface_surf/pbl_surface.f revision 309 by guez, Thu Sep 27 14:58:10 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, 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, 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    
20      ! Tout ce qui a trait aux traceurs est dans phytrac maintenant.      ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
21      ! Pour l'instant le calcul de la couche limite pour les traceurs      ! de la couche limite pour les traceurs se fait avec "cltrac" et
22      ! se fait avec cltrac et ne tient pas compte de la différentiation      ! ne tient pas compte de la diff\'erentiation des sous-fractions
23      ! des sous-fractions de sol.      ! de sol.
24    
25      ! Pour pouvoir extraire les coefficients d'échanges et le vent      use cdrag_m, only: cdrag
26      ! dans la première couche, trois champs supplémentaires ont été créés :      use clqh_m, only: clqh
27      ! zcoefh, zu1 et zv1. Pour l'instant nous avons moyenné les valeurs      use clvent_m, only: clvent
28      ! de ces trois champs sur les 4 sous-surfaces du modèle. Dans l'avenir      use coef_diff_turb_m, only: coef_diff_turb
29      ! si les informations des sous-surfaces doivent être prises en compte      USE conf_gcm_m, ONLY: lmt_pas
30      ! il faudra sortir ces mêmes champs en leur ajoutant une dimension,      USE conf_phys_m, ONLY: iflag_pbl
31      ! c'est a dire nbsrf (nombre de sous-surfaces).      USE dimphy, ONLY: klev, klon
32        USE dimsoil, ONLY: nsoilmx
     ! Auteur Z.X. Li (LMD/CNRS) date: 1993/08/18  
     ! Objet : interface de "couche limite" (diffusion verticale)  
   
     ! 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        ! 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):: play(klon, klev) ! pression au milieu de couche (Pa)
61        REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse
62        REAL, INTENT(inout):: fqsurf(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(:, :), 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
     REAL t(klon, klev), q(klon, klev)  
     REAL u(klon, klev), v(klon, klev)  
     REAL, INTENT (IN) :: paprs(klon, klev+1)  
     REAL, INTENT (IN) :: pplay(klon, klev)  
     REAL, INTENT (IN) :: rlon(klon), rlat(klon)  
     REAL cufi(klon), cvfi(klon)  
     REAL d_t(klon, klev), d_q(klon, klev)  
     REAL d_u(klon, klev), d_v(klon, klev)  
     REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)  
     REAL dflux_t(klon), dflux_q(klon)  
     !IM "slab" ocean  
     REAL flux_o(klon), flux_g(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)  
125    
126      REAL fluxlat(klon, nbsrf)      real ffonte(klon, nbsrf) ! flux thermique utilise pour fondre la neige
127        REAL, intent(inout):: run_off_lic_0(:) ! (klon)
128    
129      REAL rain_f(klon), snow_f(klon)      REAL, intent(out):: albsol(:) ! (klon)
130      REAL fder(klon)      ! albedo du sol total, visible, moyen par maille
131    
132      REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)      REAL, intent(in):: sollw(:) ! (klon)
133      REAL rugos(klon, nbsrf)      ! surface net downward longwave flux, in W m-2
     ! la nouvelle repartition des surfaces sortie de l'interface  
     REAL pctsrf_new(klon, nbsrf)  
134    
135      REAL zcoefh(klon, klev)      REAL, intent(in):: solsw(:) ! (klon)
136      REAL zu1(klon)      ! surface net downward shortwave flux, in W m-2
     REAL zv1(klon)  
137    
138      !$$$ PB ajout pour soil      REAL, intent(in):: tsol(:) ! (klon)
     LOGICAL, INTENT (IN) :: soil_model  
     !IM ajout seuils cdrm, cdrh  
     REAL cdmmax, cdhmax  
139    
140      REAL ksta, ksta_ter      ! Local:
     LOGICAL ok_kzmin  
141    
142      REAL ftsoil(klon, nsoilmx, nbsrf)      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
143      REAL ytsoil(klon, nsoilmx)      REAL fsolsw(klon, nbsrf) ! flux solaire absorb\'e pour chaque sous-surface
     REAL qsol(klon)  
144    
145      EXTERNAL clqh, clvent, coefkz, calbeta, cltrac      ! la nouvelle repartition des surfaces sortie de l'interface
146        REAL, save:: pctsrf_new_oce(klon)
147        REAL, save:: pctsrf_new_sic(klon)
148    
149      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      REAL y_fqcalving(klon), y_ffonte(klon)
150        real y_run_off_lic_0(klon), y_run_off_lic(klon)
151        REAL run_off_lic(klon) ! ruissellement total
152        REAL rugmer(klon)
153        REAL ytsoil(klon, nsoilmx)
154        REAL yts(klon), ypctsrf(klon), yz0_new(klon)
155        real yrugos(klon) ! longueur de rugosite (en m)
156      REAL yalb(klon)      REAL yalb(klon)
157      REAL yalblw(klon)      REAL snow(klon), yqsurf(klon), yagesno(klon)
158      REAL yu1(klon), yv1(klon)      real yqsol(klon) ! column-density of water in soil, in kg m-2
159      REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)      REAL yrain_fall(klon) ! liquid water mass flux (kg / m2 / s), positive down
160      REAL yrain_f(klon), ysnow_f(klon)      REAL ysnow_fall(klon) ! solid water mass flux (kg / m2 / s), positive down
161      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)  
   
162      REAL yfluxlat(klon)      REAL yfluxlat(klon)
   
163      REAL y_d_ts(klon)      REAL y_d_ts(klon)
164      REAL y_d_t(klon, klev), y_d_q(klon, klev)      REAL y_d_t(klon, klev), y_d_q(klon, klev)
165      REAL y_d_u(klon, klev), y_d_v(klon, klev)      REAL y_d_u(klon, klev), y_d_v(klon, klev)
166      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)      REAL y_flux_t(klon), y_flux_q(klon)
167      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)      REAL y_flux_u(klon), y_flux_v(klon)
168      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
169      REAL ycoefh(klon, klev), ycoefm(klon, klev)      REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
170        real ycdragh(klon), ycdragm(klon)
171      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
172      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
173      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)      REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
174        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)  
175      REAL delp(klon, klev)      REAL delp(klon, klev)
176      INTEGER i, k, nsrf      INTEGER i, k, nsrf
   
177      INTEGER ni(klon), knon, j      INTEGER ni(klon), knon, j
     ! Introduction d'une variable "pourcentage potentiel" pour tenir compte  
     ! des eventuelles apparitions et/ou disparitions de la glace de mer  
     REAL pctsrf_pot(klon, nbsrf)  
   
     REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.  
178    
179      ! maf pour sorties IOISPL en cas de debugagage      REAL pctsrf_pot(klon, nbsrf)
180        ! "pourcentage potentiel" pour tenir compte des \'eventuelles
181        ! apparitions ou disparitions de la glace de mer
182    
183      CHARACTER (80) cldebug      REAL yt2m(klon), yq2m(klon), wind10m(klon)
184      SAVE cldebug      REAL ustar(klon)
     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  
185    
186      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)  
187      REAL ypblh(klon)      REAL ypblh(klon)
188      REAL ylcl(klon)      REAL ylcl(klon)
189      REAL ycapcl(klon)      REAL ycapcl(klon)
# Line 275  contains Line 191  contains
191      REAL ycteicl(klon)      REAL ycteicl(klon)
192      REAL ypblt(klon)      REAL ypblt(klon)
193      REAL ytherm(klon)      REAL ytherm(klon)
194      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)  
195      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
196      REAL psfce(klon), patm(klon)      REAL psfce(klon), patm(klon)
197        REAL zgeo1(klon)
     REAL qairsol(klon), zgeo1(klon)  
198      REAL rugo1(klon)      REAL rugo1(klon)
199        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'  
200    
201      !------------------------------------------------------------      !------------------------------------------------------------
202    
203      ! initialisation Anne      albsol = sum(falbe * pctsrf, dim = 2)
     ytherm = 0.  
204    
205      IF (debugindex .AND. first_appel) THEN      ! R\'epartition sous maille des flux longwave et shortwave
206         first_appel = .FALSE.      ! R\'epartition du longwave par sous-surface lin\'earis\'ee
207    
208         ! initialisation sorties netcdf      forall (nsrf = 1:nbsrf)
209           fsollw(:, nsrf) = sollw + 4. * RSIGMA * tsol**3 &
210                * (tsol - ftsol(:, nsrf))
211           fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
212        END forall
213    
214         idayref = day_ini      ytherm = 0.
        CALL ymds2ju(annee_ref, 1, idayref, 0.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  
215    
216      DO k = 1, klev ! epaisseur de couche      DO k = 1, klev ! epaisseur de couche
217         DO i = 1, klon         DO i = 1, klon
218            delp(i, k) = paprs(i, k) - paprs(i, k+1)            delp(i, k) = paprs(i, k) - paprs(i, k + 1)
219         END DO         END DO
220      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  
221    
222      ! initialisation:      ! Initialization:
223        rugmer = 0.
224      DO i = 1, klon      cdragh = 0.
225         rugmer(i) = 0.0      cdragm = 0.
226         cdragh(i) = 0.0      dflux_t = 0.
227         cdragm(i) = 0.0      dflux_q = 0.
228         dflux_t(i) = 0.0      yrugos = 0.
229         dflux_q(i) = 0.0      ypaprs = 0.
230         zu1(i) = 0.0      ypplay = 0.
231         zv1(i) = 0.0      ydelp = 0.
     END DO  
     ypct = 0.0  
     yts = 0.0  
     ysnow = 0.0  
     yqsurf = 0.0  
     yalb = 0.0  
     yalblw = 0.0  
     yrain_f = 0.0  
     ysnow_f = 0.0  
     yfder = 0.0  
     ytaux = 0.0  
     ytauy = 0.0  
     ysolsw = 0.0  
     ysollw = 0.0  
     ysollwdown = 0.0  
     yrugos = 0.0  
     yu1 = 0.0  
     yv1 = 0.0  
     yrads = 0.0  
     ypaprs = 0.0  
     ypplay = 0.0  
     ydelp = 0.0  
     yu = 0.0  
     yv = 0.0  
     yt = 0.0  
     yq = 0.0  
     pctsrf_new = 0.0  
     y_flux_u = 0.0  
     y_flux_v = 0.0  
     !$$ PB  
     y_dflux_t = 0.0  
     y_dflux_q = 0.0  
     ytsoil = 999999.  
232      yrugoro = 0.      yrugoro = 0.
233      ! -- LOOP      d_ts = 0.
     yu10mx = 0.0  
     yu10my = 0.0  
     ywindsp = 0.0  
     ! -- LOOP  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           d_ts(i, nsrf) = 0.0  
        END DO  
     END DO  
     !§§§ PB  
     yfluxlat = 0.  
234      flux_t = 0.      flux_t = 0.
235      flux_q = 0.      flux_q = 0.
236      flux_u = 0.      flux_u = 0.
237      flux_v = 0.      flux_v = 0.
238      DO k = 1, klev      fluxlat = 0.
239         DO i = 1, klon      d_t = 0.
240            d_t(i, k) = 0.0      d_q = 0.
241            d_q(i, k) = 0.0      d_u = 0.
242            d_u(i, k) = 0.0      d_v = 0.
243            d_v(i, k) = 0.0      coefh = 0.
244            zcoefh(i, k) = 0.0      fqcalving = 0.
245         END DO      run_off_lic = 0.
246      END DO  
247        ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
248        ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
249        ! (\`a affiner).
250    
251      ! Boucler sur toutes les sous-fractions du sol:      pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
252        pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
     ! 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  
253      pctsrf_pot(:, is_oce) = 1. - zmasq      pctsrf_pot(:, is_oce) = 1. - zmasq
254      pctsrf_pot(:, is_sic) = 1. - zmasq      pctsrf_pot(:, is_sic) = 1. - zmasq
255    
256      DO nsrf = 1, nbsrf      ! Tester si c'est le moment de lire le fichier:
257         ! chercher les indices:      if (mod(itap - 1, lmt_pas) == 0) then
258           CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
259        endif
260    
261        ! Boucler sur toutes les sous-fractions du sol:
262    
263        loop_surface: DO nsrf = 1, nbsrf
264           ! Define ni and knon:
265    
266         ni = 0         ni = 0
267         knon = 0         knon = 0
268    
269         DO i = 1, klon         DO i = 1, klon
270            ! pour determiner le domaine a traiter on utilise les surfaces            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
271            ! "potentielles"            ! "potentielles"
272            IF (pctsrf_pot(i, nsrf) > epsfra) THEN            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
273               knon = knon + 1               knon = knon + 1
# Line 435  contains Line 275  contains
275            END IF            END IF
276         END DO         END DO
277    
278         ! variables pour avoir une sortie IOIPSL des INDEX         if_knon: IF (knon /= 0) then
279         IF (debugindex) THEN            ypctsrf(:knon) = pctsrf(ni(:knon), nsrf)
280            tabindx = 0.            yts(:knon) = ftsol(ni(:knon), nsrf)
281            DO i = 1, knon            snow(:knon) = fsnow(ni(:knon), nsrf)
282               tabindx(i) = real(i)            yqsurf(:knon) = fqsurf(ni(:knon), nsrf)
283            END DO            yalb(:knon) = falbe(ni(:knon), nsrf)
284            debugtab = 0.            yrain_fall(:knon) = rain_fall(ni(:knon))
285            ndexbg = 0            ysnow_fall(:knon) = snow_fall(ni(:knon))
286            CALL gath2cpl(tabindx, debugtab, klon, knon, iim, jjm, ni)            yagesno(:knon) = agesno(ni(:knon), nsrf)
287            CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)            yrugos(:knon) = frugs(ni(:knon), nsrf)
288         END IF            yrugoro(:knon) = rugoro(ni(:knon))
289              radsol(:knon) = fsolsw(ni(:knon), nsrf) + fsollw(ni(:knon), nsrf)
290         IF (knon==0) CYCLE            ypaprs(:knon, klev + 1) = paprs(ni(:knon), klev + 1)
291              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)  
        !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  
292    
293         !IM cf JLD : on seuille ycoefm et ycoefh            ! For continent, copy soil water content
294         IF (nsrf==is_oce) THEN            IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
           DO j = 1, knon  
              !           ycoefm(j, 1)=min(ycoefm(j, 1), 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  
295    
296         !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)  
297    
           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  
298            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  
299               DO j = 1, knon               DO j = 1, knon
300                  i = ni(j)                  i = ni(j)
301                  yq2(j, k) = q2(i, k, nsrf)                  ypaprs(j, k) = paprs(i, k)
302                    ypplay(j, k) = play(i, k)
303                    ydelp(j, k) = delp(i, k)
304                    yu(j, k) = u(i, k)
305                    yv(j, k) = v(i, k)
306                    yt(j, k) = t(i, k)
307                    yq(j, k) = q(i, k)
308               END DO               END DO
309            END DO            END DO
310    
311            !   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  
312    
313            !   iflag_pbl peut etre utilise comme longuer de melange            zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
314                   + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
315    
316            IF (iflag_pbl>=11) THEN            DO k = 2, klev
317               CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &               zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
318                    yu, yv, yteta, y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, &                    * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
319                    iflag_pbl)                    * (ypplay(:knon, k - 1) - ypplay(:knon, k))
320            ELSE            ENDDO
321               CALL yamada4(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, &  
322                    yv, yteta, y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)            CALL cdrag(nsrf, sqrt(yu(:knon, 1)**2 + yv(:knon, 1)**2), &
323                   yt(:knon, 1), yq(:knon, 1), zgeop(:knon, 1), ypaprs(:knon, 1), &
324                   yts(:knon), yqsurf(:knon), yrugos(:knon), ycdragm(:knon), &
325                   ycdragh(:knon))
326    
327              IF (iflag_pbl == 1) THEN
328                 ycdragm(:knon) = max(ycdragm(:knon), 0.)
329                 ycdragh(:knon) = max(ycdragh(:knon), 0.)
330              end IF
331    
332              ! on met un seuil pour ycdragm et ycdragh
333              IF (nsrf == is_oce) THEN
334                 ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
335                 ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
336            END IF            END IF
337    
338            ycoefm(1:knon, 1) = y_cd_m(1:knon)            IF (iflag_pbl >= 6) yq2(:knon, :) = q2(ni(:knon), :, nsrf)
339            ycoefh(1:knon, 1) = y_cd_h(1:knon)            call coef_diff_turb(nsrf, ni(:knon), ypaprs(:knon, :), &
340            ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)                 ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
341            ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)                 yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
342         END IF                 ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
343    
344         ! calculer la diffusion des vitesses "u" et "v"            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
345         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &                 ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
346              ydelp, y_d_u, y_flux_u)                 ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
347         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &                 y_flux_u(:knon))
348              ydelp, y_d_v, y_flux_v)            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
349                   ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
350         ! pour le couplage                 ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
351         ytaux = y_flux_u(:, 1)                 y_flux_v(:knon))
352         ytauy = y_flux_v(:, 1)  
353              CALL clqh(julien, nsrf, ni(:knon), ytsoil(:knon, :), yqsol(:knon), &
354         ! FH modif sur le cdrag temperature                 mu0(ni(:knon)), yrugos(:knon), yrugoro(:knon), yu(:knon, 1), &
355         !$$$PB : déplace dans clcdrag                 yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), yt(:knon, :), &
356         !$$$      do i=1, knon                 yq(:knon, :), yts(:knon), ypaprs(:knon, :), ypplay(:knon, :), &
357         !$$$         ycoefh(i, 1)=ycoefm(i, 1)*0.8                 ydelp(:knon, :), radsol(:knon), yalb(:knon), snow(:knon), &
358         !$$$      enddo                 yqsurf(:knon), yrain_fall(:knon), ysnow_fall(:knon), &
359                   yfluxlat(:knon), pctsrf_new_sic(ni(:knon)), yagesno(:knon), &
360         ! calculer la diffusion de "q" et de "h"                 y_d_t(:knon, :), y_d_q(:knon, :), y_d_ts(:knon), &
361         CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,&                 yz0_new(:knon), y_flux_t(:knon), y_flux_q(:knon), &
362              cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil,&                 y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving(:knon), &
363              yqsol, ok_veget, ocean, npas, nexca, rmu0, co2_ppm, yrugos,&                 y_ffonte(:knon), y_run_off_lic_0(:knon), y_run_off_lic(:knon))
             yrugoro, yu1, yv1, ycoefh, yt, yq, yts, ypaprs, ypplay,&  
             ydelp, yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, &  
             yfder, ytaux, ytauy, ywindsp, ysollw, ysollwdown, ysolsw,&  
             yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts,&  
             yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,&  
             y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g,&  
             ytslab, y_seaice)  
   
        ! calculer la longueur de rugosite sur ocean  
        yrugm = 0.  
        IF (nsrf==is_oce) THEN  
           DO j = 1, knon  
              yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &  
                   0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2))  
              yrugm(j) = max(1.5E-05, yrugm(j))  
           END DO  
        END IF  
        DO j = 1, knon  
           y_dflux_t(j) = y_dflux_t(j)*ypct(j)  
           y_dflux_q(j) = y_dflux_q(j)*ypct(j)  
           yu1(j) = yu1(j)*ypct(j)  
           yv1(j) = yv1(j)*ypct(j)  
        END DO  
364    
365         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  
366    
367         evap(:, nsrf) = -flux_q(:, 1, nsrf)            yrugm = 0.
368    
369         albe(:, nsrf) = 0.            IF (nsrf == is_oce) THEN
370         alblw(:, nsrf) = 0.               DO j = 1, knon
371         snow(:, nsrf) = 0.                  yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
372         qsurf(:, nsrf) = 0.                       / rg + 0.11 * 14E-6 &
373         rugos(:, nsrf) = 0.                       / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
374         fluxlat(:, nsrf) = 0.                  yrugm(j) = max(1.5E-05, yrugm(j))
375         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)  
376            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  
377    
        DO j = 1, knon  
           i = ni(j)  
378            DO k = 1, klev            DO k = 1, klev
379               d_t(i, k) = d_t(i, k) + y_d_t(j, k)               DO j = 1, knon
380               d_q(i, k) = d_q(i, k) + y_d_q(j, k)                  i = ni(j)
381               !$$$ PB        flux_t(i, k) = flux_t(i, k) + y_flux_t(j, k)                  y_d_t(j, k) = y_d_t(j, k) * ypctsrf(j)
382               !$$$         flux_q(i, k) = flux_q(i, k) + y_flux_q(j, k)                  y_d_q(j, k) = y_d_q(j, k) * ypctsrf(j)
383               d_u(i, k) = d_u(i, k) + y_d_u(j, k)                  y_d_u(j, k) = y_d_u(j, k) * ypctsrf(j)
384               d_v(i, k) = d_v(i, k) + y_d_v(j, k)                  y_d_v(j, k) = y_d_v(j, k) * ypctsrf(j)
385               !$$$  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)  
386            END DO            END DO
        END DO  
   
        !cc diagnostic t, q a 2m et u, v a 10m  
387    
388         DO j = 1, knon            flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
389            i = ni(j)            flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
390            uzon(j) = yu(j, 1) + y_d_u(j, 1)            flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
391            vmer(j) = yv(j, 1) + y_d_v(j, 1)            flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
392            tair1(j) = yt(j, 1) + y_d_t(j, 1)  
393            qair1(j) = yq(j, 1) + y_d_q(j, 1)            falbe(:, nsrf) = 0.
394            zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &            fsnow(:, nsrf) = 0.
395                 1)))*(ypaprs(j, 1)-ypplay(j, 1))            fqsurf(:, nsrf) = 0.
396            tairsol(j) = yts(j) + y_d_ts(j)            frugs(:, nsrf) = 0.
397            rugo1(j) = yrugos(j)            DO j = 1, knon
398            IF (nsrf==is_oce) THEN               i = ni(j)
399               rugo1(j) = rugos(i, nsrf)               d_ts(i, nsrf) = y_d_ts(j)
400                 falbe(i, nsrf) = yalb(j)
401                 fsnow(i, nsrf) = snow(j)
402                 fqsurf(i, nsrf) = yqsurf(j)
403                 frugs(i, nsrf) = yz0_new(j)
404                 fluxlat(i, nsrf) = yfluxlat(j)
405                 IF (nsrf == is_oce) THEN
406                    rugmer(i) = yrugm(j)
407                    frugs(i, nsrf) = yrugm(j)
408                 END IF
409                 agesno(i, nsrf) = yagesno(j)
410                 fqcalving(i, nsrf) = y_fqcalving(j)
411                 ffonte(i, nsrf) = y_ffonte(j)
412                 cdragh(i) = cdragh(i) + ycdragh(j) * ypctsrf(j)
413                 cdragm(i) = cdragm(i) + ycdragm(j) * ypctsrf(j)
414                 dflux_t(i) = dflux_t(i) + y_dflux_t(j) * ypctsrf(j)
415                 dflux_q(i) = dflux_q(i) + y_dflux_q(j) * ypctsrf(j)
416              END DO
417              IF (nsrf == is_ter) THEN
418                 qsol(ni(:knon)) = yqsol(:knon)
419              else IF (nsrf == is_lic) THEN
420                 DO j = 1, knon
421                    i = ni(j)
422                    run_off_lic_0(i) = y_run_off_lic_0(j)
423                    run_off_lic(i) = y_run_off_lic(j)
424                 END DO
425            END IF            END IF
           psfce(j) = ypaprs(j, 1)  
           patm(j) = ypplay(j, 1)  
426    
427            qairsol(j) = yqsurf(j)            ftsoil(:, :, nsrf) = 0.
428         END DO            ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)
429    
430         CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &            DO j = 1, knon
431              tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &               i = ni(j)
432              yu10m, yustar)               DO k = 1, klev
433         !IM 081204 END                  d_t(i, k) = d_t(i, k) + y_d_t(j, k)
434                    d_q(i, k) = d_q(i, k) + y_d_q(j, k)
435         DO j = 1, knon                  d_u(i, k) = d_u(i, k) + y_d_u(j, k)
436            i = ni(j)                  d_v(i, k) = d_v(i, k) + y_d_v(j, k)
437            t2m(i, nsrf) = yt2m(j)               END DO
438            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)  
439    
440         END DO            forall (k = 2:klev) coefh(ni(:knon), k) &
441                   = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypctsrf(:knon)
442    
443         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  
444    
445         DO j = 1, knon            DO j = 1, knon
           DO k = 1, klev + 1  
446               i = ni(j)               i = ni(j)
447               q2(i, k, nsrf) = yq2(j, k)               u1(j) = yu(j, 1) + y_d_u(j, 1)
448                 v1(j) = yv(j, 1) + y_d_v(j, 1)
449                 tair1(j) = yt(j, 1) + y_d_t(j, 1)
450                 qair1(j) = yq(j, 1) + y_d_q(j, 1)
451                 zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
452                      1))) * (ypaprs(j, 1)-ypplay(j, 1))
453                 tairsol(j) = yts(j) + y_d_ts(j)
454                 rugo1(j) = yrugos(j)
455                 IF (nsrf == is_oce) THEN
456                    rugo1(j) = frugs(i, nsrf)
457                 END IF
458                 psfce(j) = ypaprs(j, 1)
459                 patm(j) = ypplay(j, 1)
460            END DO            END DO
461         END DO  
462         !IM "slab" ocean            CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), qair1, &
463         IF (nsrf==is_oce) THEN                 zgeo1, tairsol, yqsurf(:knon), rugo1, psfce, patm, yt2m, yq2m, &
464                   yt10m, yq10m, wind10m(:knon), ustar(:knon))
465    
466            DO j = 1, knon            DO j = 1, knon
              ! on projette sur la grille globale  
467               i = ni(j)               i = ni(j)
468               IF (pctsrf_new(i, is_oce)>epsfra) THEN               t2m(i, nsrf) = yt2m(j)
469                  flux_o(i) = y_flux_o(j)               q2m(i, nsrf) = yq2m(j)
470               ELSE  
471                  flux_o(i) = 0.               u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
472               END IF                    / sqrt(u1(j)**2 + v1(j)**2)
473                 v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
474                      / sqrt(u1(j)**2 + v1(j)**2)
475            END DO            END DO
        END IF  
476    
477         IF (nsrf==is_sic) THEN            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
478                   y_flux_q(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
479                   yq(:knon, :), ypblh(:knon), ycapcl, yoliqcl, ycteicl, ypblt, &
480                   ytherm, ylcl)
481    
482            DO j = 1, knon            DO j = 1, knon
483               i = ni(j)               i = ni(j)
484               ! On pondère lorsque l'on fait le bilan au sol :               pblh(i, nsrf) = ypblh(j)
485               ! flux_g(i) = y_flux_g(j)*ypct(j)               plcl(i, nsrf) = ylcl(j)
486               IF (pctsrf_new(i, is_sic)>epsfra) THEN               capcl(i, nsrf) = ycapcl(j)
487                  flux_g(i) = y_flux_g(j)               oliqcl(i, nsrf) = yoliqcl(j)
488               ELSE               cteicl(i, nsrf) = ycteicl(j)
489                  flux_g(i) = 0.               pblt(i, nsrf) = ypblt(j)
490               END IF               therm(i, nsrf) = ytherm(j)
491            END DO            END DO
492    
493         END IF            IF (iflag_pbl >= 6) q2(ni(:knon), :, nsrf) = yq2(:knon, :)
494         !nsrf.EQ.is_sic                                                     else
495         IF (ocean=='slab  ') THEN            fsnow(:, nsrf) = 0.
496            IF (nsrf==is_oce) THEN         end IF if_knon
497               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  
498    
499      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
500      ! A rajouter: conservation de l'albedo      frugs(:, is_oce) = rugmer
501        pctsrf(:, is_oce) = pctsrf_new_oce
502        pctsrf(:, is_sic) = pctsrf_new_sic
503    
504      rugos(:, is_oce) = rugmer      CALL histwrite_phy("run_off_lic", run_off_lic)
     pctsrf = pctsrf_new  
505    
506    END SUBROUTINE clmain    END SUBROUTINE pbl_surface
507    
508  end module clmain_m  end module pbl_surface_m
Legend:
Removed from v.38  
changed lines
  Added in v.309
  ViewVC Help
Powered by ViewVC 1.1.21