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

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

  ViewVC Help
Powered by ViewVC 1.1.21