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

Diff of /trunk/phylmd/pbl_surface.f

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

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