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