/[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 72 by guez, Tue Jul 23 13:00:07 2013 UTC trunk/phylmd/pbl_surface.f revision 276 by guez, Thu Jul 12 14:49:20 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, 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, 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, fder, rlon, rlat, &         q2, dflux_t, dflux_q, coefh, t2m, q2m, u10m_srf, v10m_srf, pblh, capcl, &
12         rugos, debut, 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 cdrag_m, only: cdrag
     ! 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
29      USE dimens_m, ONLY: iim, jjm      USE dimphy, ONLY: klev, klon
     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 phyetat0_m, only: zmasq
35      use ustarhb_m, only: ustarhb      use stdlevvar_m, only: stdlevvar
36      use vdif_kcay_m, only: vdif_kcay      USE suphec_m, ONLY: rd, rg
37      use yamada4_m, only: yamada4      use time_phylmdz, only: itap
   
     ! Arguments:  
38    
39      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(inout):: pctsrf(klon, nbsrf)  
40    
41      ! la nouvelle repartition des surfaces sortie de l'interface      REAL, INTENT(inout):: pctsrf(klon, nbsrf)
42      REAL, INTENT(out):: pctsrf_new(klon, nbsrf)      ! tableau des pourcentages de surface de chaque maille
43    
44      REAL, INTENT(IN):: t(klon, klev) ! temperature (K)      REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
45      REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg/kg)      REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg / kg)
46      REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse      REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
47      INTEGER, INTENT(IN):: jour ! jour de l'annee en cours      INTEGER, INTENT(IN):: julien ! jour de l'annee en cours
48      REAL, intent(in):: rmu0(klon) ! cosinus de l'angle solaire zenithal          REAL, intent(in):: mu0(klon) ! cosinus de l'angle solaire zenithal    
49      REAL co2_ppm ! taux CO2 atmosphere      REAL, INTENT(IN):: ftsol(:, :) ! (klon, nbsrf) temp\'erature du sol (en K)
     LOGICAL ok_veget  
     CHARACTER(len=*), INTENT(IN):: ocean  
     REAL ts(klon, nbsrf) ! input-R- temperature du sol (en Kelvin)  
     LOGICAL, INTENT(IN):: soil_model  
50      REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh      REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
51      REAL ksta, ksta_ter  
52      LOGICAL ok_kzmin      REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
53      REAL ftsoil(klon, nsoilmx, nbsrf)      ! soil temperature of surface fraction
54      REAL qsol(klon)  
55      REAL, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)      REAL, INTENT(inout):: qsol(:) ! (klon)
56        ! column-density of water in soil, in kg m-2
57    
58        REAL, INTENT(IN):: paprs(klon, klev + 1) ! pression a intercouche (Pa)
59      REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)      REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
60      REAL snow(klon, nbsrf)      REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse
61      REAL qsurf(klon, nbsrf)      REAL qsurf(klon, nbsrf)
62      REAL evap(klon, nbsrf)      REAL evap(klon, nbsrf)
63      REAL albe(klon, nbsrf)      REAL, intent(inout):: falbe(klon, nbsrf)
64      REAL alblw(klon, nbsrf)      REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)
   
     REAL fluxlat(klon, nbsrf)  
65    
66      REAL, intent(in):: rain_fall(klon), snow_f(klon)      REAL, intent(in):: rain_fall(klon)
67      REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)      ! 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  
68    
69      REAL rugos(klon, nbsrf)      REAL, intent(in):: snow_f(klon)
70      ! rugos----input-R- longeur de rugosite (en m)      ! solid water mass flux (kg / m2 / s), positive down
71    
72      LOGICAL, INTENT(IN):: debut      REAL, INTENT(IN):: fsolsw(klon, nbsrf), fsollw(klon, nbsrf)
73        REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)
74      real agesno(klon, nbsrf)      real agesno(klon, nbsrf)
75      REAL, INTENT(IN):: rugoro(klon)      REAL, INTENT(IN):: rugoro(klon)
76    
# Line 108  contains Line 81  contains
81      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
82      ! changement pour "u" et "v"      ! changement pour "u" et "v"
83    
84      REAL d_ts(klon, nbsrf)      REAL, intent(out):: d_ts(:, :) ! (klon, nbsrf) variation of ftsol
     ! d_ts-----output-R- le changement pour "ts"  
85    
86      REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)      REAL, intent(out):: flux_t(klon, nbsrf)
87      ! 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
88      !                    (orientation positive vers le bas)      ! le bas) à la surface
89      ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)  
90        REAL, intent(out):: flux_q(klon, nbsrf)
91      REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)      ! flux de vapeur d'eau (kg / m2 / s) à la surface
92      ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal  
93      ! 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)
94        ! tension du vent (flux turbulent de vent) à la surface, en Pa
95    
96      REAL, INTENT(out):: cdragh(klon), cdragm(klon)      REAL, INTENT(out):: cdragh(klon), cdragm(klon)
97      real q2(klon, klev+1, nbsrf)      real q2(klon, klev + 1, nbsrf)
98    
99      REAL dflux_t(klon), dflux_q(klon)      REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)
100      ! dflux_t derive du flux sensible      ! dflux_t derive du flux sensible
101      ! dflux_q derive du flux latent      ! dflux_q derive du flux latent
102      !IM "slab" ocean      ! IM "slab" ocean
103    
104      REAL, intent(out):: ycoefh(klon, klev)      REAL, intent(out):: coefh(:, 2:) ! (klon, 2:klev)
105      REAL, intent(out):: zu1(klon)      ! Pour pouvoir extraire les coefficients d'\'echange, le champ
106      REAL zv1(klon)      ! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
107      REAL t2m(klon, nbsrf), q2m(klon, nbsrf)      ! ce champ sur les quatre sous-surfaces du mod\`ele.
108      REAL u10m(klon, nbsrf), v10m(klon, nbsrf)  
109        REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)
110      !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds  
111      ! physiq ce qui permet de sortir les grdeurs par sous surface)      REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)
112      REAL pblh(klon, nbsrf)      ! composantes du vent \`a 10m sans spirale d'Ekman
113      ! pblh------- HCL  
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)      REAL capcl(klon, nbsrf)
119      REAL oliqcl(klon, nbsrf)      REAL oliqcl(klon, nbsrf)
120      REAL cteicl(klon, nbsrf)      REAL cteicl(klon, nbsrf)
121      REAL pblt(klon, nbsrf)      REAL, INTENT(inout):: pblt(klon, nbsrf) ! T au nveau HCL
     ! pblT------- T au nveau HCL  
122      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  
123      REAL plcl(klon, nbsrf)      REAL plcl(klon, nbsrf)
124      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
125      ! ffonte----Flux thermique utilise pour fondre la neige      ! ffonte----Flux thermique utilise pour fondre la neige
126      ! 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
127      !           hauteur de neige, en kg/m2/s      !           hauteur de neige, en kg / m2 / s
128      REAL run_off_lic_0(klon)      REAL run_off_lic_0(klon)
129    
130      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  
131    
132      REAL seaice(klon)      LOGICAL:: firstcal = .true.
     ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')  
133    
134      ! Local:      ! la nouvelle repartition des surfaces sortie de l'interface
135        REAL, save:: pctsrf_new_oce(klon)
136        REAL, save:: pctsrf_new_sic(klon)
137    
     REAL y_flux_o(klon), y_flux_g(klon)  
     real ytslab(klon)  
     real y_seaice(klon)  
138      REAL y_fqcalving(klon), y_ffonte(klon)      REAL y_fqcalving(klon), y_ffonte(klon)
139      real y_run_off_lic_0(klon)      real y_run_off_lic_0(klon)
   
140      REAL rugmer(klon)      REAL rugmer(klon)
   
141      REAL ytsoil(klon, nsoilmx)      REAL ytsoil(klon, nsoilmx)
142        REAL yts(klon), ypct(klon), yz0_new(klon)
143      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      real yrugos(klon) ! longeur de rugosite (en m)
144      REAL yalb(klon)      REAL yalb(klon)
145      REAL yalblw(klon)      REAL snow(klon), yqsurf(klon), yagesno(klon)
146      REAL yu1(klon), yv1(klon)      real yqsol(klon) ! column-density of water in soil, in kg m-2
147      ! 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
148      ! 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)  
     REAL yfder(klon), ytaux(klon), ytauy(klon)  
149      REAL yrugm(klon), yrads(klon), yrugoro(klon)      REAL yrugm(klon), yrads(klon), yrugoro(klon)
   
150      REAL yfluxlat(klon)      REAL yfluxlat(klon)
   
151      REAL y_d_ts(klon)      REAL y_d_ts(klon)
152      REAL y_d_t(klon, klev), y_d_q(klon, klev)      REAL y_d_t(klon, klev), y_d_q(klon, klev)
153      REAL y_d_u(klon, klev), y_d_v(klon, klev)      REAL y_d_u(klon, klev), y_d_v(klon, klev)
154      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)      REAL y_flux_t(klon), y_flux_q(klon)
155      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)      REAL y_flux_u(klon), y_flux_v(klon)
156      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
157      REAL coefh(klon, klev), coefm(klon, klev)      REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
158        real ycdragh(klon), ycdragm(klon)
159      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
160      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
161      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)      REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
162        REAL yq2(klon, klev + 1)
     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)  
163      REAL delp(klon, klev)      REAL delp(klon, klev)
164      INTEGER i, k, nsrf      INTEGER i, k, nsrf
   
165      INTEGER ni(klon), knon, j      INTEGER ni(klon), knon, j
166    
167      REAL pctsrf_pot(klon, nbsrf)      REAL pctsrf_pot(klon, nbsrf)
168      ! "pourcentage potentiel" pour tenir compte des éventuelles      ! "pourcentage potentiel" pour tenir compte des \'eventuelles
169      ! apparitions ou disparitions de la glace de mer      ! apparitions ou disparitions de la glace de mer
170    
171      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.      REAL yt2m(klon), yq2m(klon), wind10m(klon)
172        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  
173    
174      REAL yt10m(klon), yq10m(klon)      REAL yt10m(klon), yq10m(klon)
175      REAL ypblh(klon)      REAL ypblh(klon)
# Line 260  contains Line 179  contains
179      REAL ycteicl(klon)      REAL ycteicl(klon)
180      REAL ypblt(klon)      REAL ypblt(klon)
181      REAL ytherm(klon)      REAL ytherm(klon)
182      REAL ytrmb1(klon)      REAL u1(klon), v1(klon)
     REAL ytrmb2(klon)  
     REAL ytrmb3(klon)  
     REAL uzon(klon), vmer(klon)  
183      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
184      REAL psfce(klon), patm(klon)      REAL psfce(klon), patm(klon)
185    
186      REAL qairsol(klon), zgeo1(klon)      REAL qairsol(klon), zgeo1(klon)
187      REAL rugo1(klon)      REAL rugo1(klon)
188        REAL zgeop(klon, klev)
     ! utiliser un jeu de fonctions simples                
     LOGICAL zxli  
     PARAMETER (zxli=.FALSE.)  
189    
190      !------------------------------------------------------------      !------------------------------------------------------------
191    
192      ytherm = 0.      ytherm = 0.
193    
     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  
   
194      DO k = 1, klev ! epaisseur de couche      DO k = 1, klev ! epaisseur de couche
195         DO i = 1, klon         DO i = 1, klon
196            delp(i, k) = paprs(i, k) - paprs(i, k+1)            delp(i, k) = paprs(i, k) - paprs(i, k + 1)
197         END DO         END DO
198      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  
199    
200      ! Initialization:      ! Initialization:
201      rugmer = 0.      rugmer = 0.
# Line 325  contains Line 203  contains
203      cdragm = 0.      cdragm = 0.
204      dflux_t = 0.      dflux_t = 0.
205      dflux_q = 0.      dflux_q = 0.
     zu1 = 0.  
     zv1 = 0.  
206      ypct = 0.      ypct = 0.
     yts = 0.  
     ysnow = 0.  
207      yqsurf = 0.      yqsurf = 0.
     yalb = 0.  
     yalblw = 0.  
208      yrain_f = 0.      yrain_f = 0.
209      ysnow_f = 0.      ysnow_f = 0.
     yfder = 0.  
     ytaux = 0.  
     ytauy = 0.  
     ysolsw = 0.  
     ysollw = 0.  
210      yrugos = 0.      yrugos = 0.
     yu1 = 0.  
     yv1 = 0.  
     yrads = 0.  
211      ypaprs = 0.      ypaprs = 0.
212      ypplay = 0.      ypplay = 0.
213      ydelp = 0.      ydelp = 0.
# Line 351  contains Line 215  contains
215      yv = 0.      yv = 0.
216      yt = 0.      yt = 0.
217      yq = 0.      yq = 0.
     pctsrf_new = 0.  
     y_flux_u = 0.  
     y_flux_v = 0.  
     !$$ PB  
218      y_dflux_t = 0.      y_dflux_t = 0.
219      y_dflux_q = 0.      y_dflux_q = 0.
     ytsoil = 999999.  
220      yrugoro = 0.      yrugoro = 0.
     ! -- LOOP  
     yu10mx = 0.  
     yu10my = 0.  
     ywindsp = 0.  
     ! -- LOOP  
221      d_ts = 0.      d_ts = 0.
     !§§§ PB  
     yfluxlat = 0.  
222      flux_t = 0.      flux_t = 0.
223      flux_q = 0.      flux_q = 0.
224      flux_u = 0.      flux_u = 0.
225      flux_v = 0.      flux_v = 0.
226        fluxlat = 0.
227      d_t = 0.      d_t = 0.
228      d_q = 0.      d_q = 0.
229      d_u = 0.      d_u = 0.
230      d_v = 0.      d_v = 0.
231      ycoefh = 0.      coefh = 0.
   
     ! Boucler sur toutes les sous-fractions du sol:  
232    
233      ! Initialisation des "pourcentages potentiels". On considère ici qu'on      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
234      ! peut avoir potentiellement de la glace sur tout le domaine océanique      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
235      ! (à affiner)      ! (\`a affiner)
236    
237      pctsrf_pot = pctsrf      pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
238        pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
239      pctsrf_pot(:, is_oce) = 1. - zmasq      pctsrf_pot(:, is_oce) = 1. - zmasq
240      pctsrf_pot(:, is_sic) = 1. - zmasq      pctsrf_pot(:, is_sic) = 1. - zmasq
241    
242        ! Tester si c'est le moment de lire le fichier:
243        if (mod(itap - 1, lmt_pas) == 0) then
244           CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
245        endif
246    
247        ! Boucler sur toutes les sous-fractions du sol:
248    
249      loop_surface: DO nsrf = 1, nbsrf      loop_surface: DO nsrf = 1, nbsrf
250         ! Chercher les indices :         ! Chercher les indices :
251         ni = 0         ni = 0
252         knon = 0         knon = 0
253         DO i = 1, klon         DO i = 1, klon
254            ! Pour déterminer le domaine à traiter, on utilise les surfaces            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
255            ! "potentielles"            ! "potentielles"
256            IF (pctsrf_pot(i, nsrf) > epsfra) THEN            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
257               knon = knon + 1               knon = knon + 1
# Line 400  contains Line 259  contains
259            END IF            END IF
260         END DO         END DO
261    
        ! 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  
   
262         if_knon: IF (knon /= 0) then         if_knon: IF (knon /= 0) then
263            DO j = 1, knon            DO j = 1, knon
264               i = ni(j)               i = ni(j)
265               ypct(j) = pctsrf(i, nsrf)               ypct(j) = pctsrf(i, nsrf)
266               yts(j) = ts(i, nsrf)               yts(j) = ftsol(i, nsrf)
267               ytslab(i) = tslab(i)               snow(j) = fsnow(i, nsrf)
              ysnow(j) = snow(i, nsrf)  
268               yqsurf(j) = qsurf(i, nsrf)               yqsurf(j) = qsurf(i, nsrf)
269               yalb(j) = albe(i, nsrf)               yalb(j) = falbe(i, nsrf)
              yalblw(j) = alblw(i, nsrf)  
270               yrain_f(j) = rain_fall(i)               yrain_f(j) = rain_fall(i)
271               ysnow_f(j) = snow_f(i)               ysnow_f(j) = snow_f(i)
272               yagesno(j) = agesno(i, nsrf)               yagesno(j) = agesno(i, nsrf)
273               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)  
              yrugos(j) = rugos(i, nsrf)  
274               yrugoro(j) = rugoro(i)               yrugoro(j) = rugoro(i)
275               yu1(j) = u1lay(i)               yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)
276               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)  
277               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))  
278            END DO            END DO
279    
280            ! IF bucket model for continent, copy soil water content            ! For continent, copy soil water content
281            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  
282    
283            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  
284    
285            DO k = 1, klev            DO k = 1, klev
286               DO j = 1, knon               DO j = 1, knon
# Line 472  contains Line 295  contains
295               END DO               END DO
296            END DO            END DO
297    
298            ! calculer Cdrag et les coefficients d'echange            ! Calculer les géopotentiels de chaque couche:
299            CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &  
300                 yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))            zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
301                   + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
302    
303              DO k = 2, klev
304                 zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
305                      * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
306                      * (ypplay(:knon, k - 1) - ypplay(:knon, k))
307              ENDDO
308    
309              CALL cdrag(nsrf, sqrt(yu(:knon, 1)**2 + yv(:knon, 1)**2), &
310                   yt(:knon, 1), yq(:knon, 1), zgeop(:knon, 1), ypaprs(:knon, 1), &
311                   yts(:knon), yqsurf(:knon), yrugos(:knon), ycdragm(:knon), &
312                   ycdragh(:knon))
313    
314            IF (iflag_pbl == 1) THEN            IF (iflag_pbl == 1) THEN
315               CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)               ycdragm(:knon) = max(ycdragm(:knon), 0.)
316               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))               ycdragh(:knon) = max(ycdragh(:knon), 0.)
317               coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))            end IF
           END IF  
318    
319            ! on met un seuil pour coefm et coefh            ! on met un seuil pour ycdragm et ycdragh
320            IF (nsrf == is_oce) THEN            IF (nsrf == is_oce) THEN
321               coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)               ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
322               coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)               ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
323            END IF            END IF
324    
325            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  
326               DO k = 1, klev + 1               DO k = 1, klev + 1
327                  DO j = 1, knon                  DO j = 1, knon
328                     i = ni(j)                     i = ni(j)
329                     yq2(j, k) = q2(i, k, nsrf)                     yq2(j, k) = q2(i, k, nsrf)
330                  END DO                  END DO
331               END DO               END DO
332              end IF
333    
334               CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)            call coef_diff_turb(dtime, nsrf, ni(:knon), ypaprs(:knon, :), &
335                   ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
336               IF (prt_level > 9) THEN                 yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
337                  PRINT *, 'USTAR = ', yustar                 ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
338               END IF  
339              CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
340               ! iflag_pbl peut être utilisé comme longueur de mélange                 ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
341                   ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
342               IF (iflag_pbl >= 11) THEN                 y_flux_u(:knon))
343                  CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &            CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
344                       yu, yv, yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, &                 ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
345                       yustar, iflag_pbl)                 ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
346               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)  
347    
348            ! calculer la diffusion de "q" et de "h"            ! calculer la diffusion de "q" et de "h"
349            CALL clqh(dtime, itap, jour, debut, rlat, knon, nsrf, ni, pctsrf, &            CALL clqh(dtime, julien, firstcal, nsrf, ni(:knon), &
350                 soil_model, ytsoil, yqsol, ok_veget, ocean, rmu0, co2_ppm, &                 ytsoil(:knon, :), yqsol(:knon), mu0, yrugos, yrugoro, &
351                 yrugos, yrugoro, yu1, yv1, coefh(:knon, :), yt, yq, yts, &                 yu(:knon, 1), yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), &
352                 ypaprs, ypplay, ydelp, yrads, yalb, yalblw, ysnow, yqsurf, &                 yt, yq, yts(:knon), ypaprs, ypplay, ydelp, yrads(:knon), &
353                 yrain_f, ysnow_f, yfder, ysolsw, yfluxlat, pctsrf_new, &                 yalb(:knon), snow(:knon), yqsurf, yrain_f, ysnow_f, &
354                 yagesno, y_d_t, y_d_q, y_d_ts, yz0_new, y_flux_t, y_flux_q, &                 yfluxlat(:knon), pctsrf_new_sic, yagesno(:knon), y_d_t, y_d_q, &
355                 y_dflux_t, y_dflux_q, y_fqcalving, y_ffonte, y_run_off_lic_0, &                 y_d_ts(:knon), yz0_new, y_flux_t(:knon), y_flux_q(:knon), &
356                 y_flux_o, y_flux_g, ytslab, y_seaice)                 y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving, y_ffonte, &
357                   y_run_off_lic_0)
358    
359            ! calculer la longueur de rugosite sur ocean            ! calculer la longueur de rugosite sur ocean
360            yrugm = 0.            yrugm = 0.
361            IF (nsrf == is_oce) THEN            IF (nsrf == is_oce) THEN
362               DO j = 1, knon               DO j = 1, knon
363                  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) &
364                       0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))                       / rg + 0.11 * 14E-6 &
365                         / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
366                  yrugm(j) = max(1.5E-05, yrugm(j))                  yrugm(j) = max(1.5E-05, yrugm(j))
367               END DO               END DO
368            END IF            END IF
369            DO j = 1, knon            DO j = 1, knon
370               y_dflux_t(j) = y_dflux_t(j)*ypct(j)               y_dflux_t(j) = y_dflux_t(j) * ypct(j)
371               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)  
372            END DO            END DO
373    
374            DO k = 1, klev            DO k = 1, klev
375               DO j = 1, knon               DO j = 1, knon
376                  i = ni(j)                  i = ni(j)
377                  coefh(j, k) = coefh(j, k)*ypct(j)                  y_d_t(j, k) = y_d_t(j, k) * ypct(j)
378                  coefm(j, k) = coefm(j, k)*ypct(j)                  y_d_q(j, k) = y_d_q(j, k) * ypct(j)
379                  y_d_t(j, k) = y_d_t(j, k)*ypct(j)                  y_d_u(j, k) = y_d_u(j, k) * ypct(j)
380                  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)  
381               END DO               END DO
382            END DO            END DO
383    
384            evap(:, nsrf) = -flux_q(:, 1, nsrf)            flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
385              flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
386              flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
387              flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
388    
389              evap(:, nsrf) = -flux_q(:, nsrf)
390    
391            albe(:, nsrf) = 0.            falbe(:, nsrf) = 0.
392            alblw(:, nsrf) = 0.            fsnow(:, nsrf) = 0.
           snow(:, nsrf) = 0.  
393            qsurf(:, nsrf) = 0.            qsurf(:, nsrf) = 0.
394            rugos(:, nsrf) = 0.            frugs(:, nsrf) = 0.
           fluxlat(:, nsrf) = 0.  
395            DO j = 1, knon            DO j = 1, knon
396               i = ni(j)               i = ni(j)
397               d_ts(i, nsrf) = y_d_ts(j)               d_ts(i, nsrf) = y_d_ts(j)
398               albe(i, nsrf) = yalb(j)               falbe(i, nsrf) = yalb(j)
399               alblw(i, nsrf) = yalblw(j)               fsnow(i, nsrf) = snow(j)
              snow(i, nsrf) = ysnow(j)  
400               qsurf(i, nsrf) = yqsurf(j)               qsurf(i, nsrf) = yqsurf(j)
401               rugos(i, nsrf) = yz0_new(j)               frugs(i, nsrf) = yz0_new(j)
402               fluxlat(i, nsrf) = yfluxlat(j)               fluxlat(i, nsrf) = yfluxlat(j)
403               IF (nsrf == is_oce) THEN               IF (nsrf == is_oce) THEN
404                  rugmer(i) = yrugm(j)                  rugmer(i) = yrugm(j)
405                  rugos(i, nsrf) = yrugm(j)                  frugs(i, nsrf) = yrugm(j)
406               END IF               END IF
407               agesno(i, nsrf) = yagesno(j)               agesno(i, nsrf) = yagesno(j)
408               fqcalving(i, nsrf) = y_fqcalving(j)               fqcalving(i, nsrf) = y_fqcalving(j)
409               ffonte(i, nsrf) = y_ffonte(j)               ffonte(i, nsrf) = y_ffonte(j)
410               cdragh(i) = cdragh(i) + coefh(j, 1)               cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
411               cdragm(i) = cdragm(i) + coefm(j, 1)               cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
412               dflux_t(i) = dflux_t(i) + y_dflux_t(j)               dflux_t(i) = dflux_t(i) + y_dflux_t(j)
413               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)  
414            END DO            END DO
415            IF (nsrf == is_ter) THEN            IF (nsrf == is_ter) THEN
416               DO j = 1, knon               qsol(ni(:knon)) = yqsol(:knon)
417                  i = ni(j)            else IF (nsrf == is_lic) THEN
                 qsol(i) = yqsol(j)  
              END DO  
           END IF  
           IF (nsrf == is_lic) THEN  
418               DO j = 1, knon               DO j = 1, knon
419                  i = ni(j)                  i = ni(j)
420                  run_off_lic_0(i) = y_run_off_lic_0(j)                  run_off_lic_0(i) = y_run_off_lic_0(j)
421               END DO               END DO
422            END IF            END IF
423            !$$$ PB ajout pour soil  
424            ftsoil(:, :, nsrf) = 0.            ftsoil(:, :, nsrf) = 0.
425            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  
426    
427            DO j = 1, knon            DO j = 1, knon
428               i = ni(j)               i = ni(j)
# Line 656  contains Line 431  contains
431                  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)
432                  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)
433                  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)  
434               END DO               END DO
435            END DO            END DO
436    
437            !cc diagnostic t, q a 2m et u, v a 10m            forall (k = 2:klev) coefh(ni(:knon), k) &
438                   = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypct(:knon)
439    
440              ! diagnostic t, q a 2m et u, v a 10m
441    
442            DO j = 1, knon            DO j = 1, knon
443               i = ni(j)               i = ni(j)
444               uzon(j) = yu(j, 1) + y_d_u(j, 1)               u1(j) = yu(j, 1) + y_d_u(j, 1)
445               vmer(j) = yv(j, 1) + y_d_v(j, 1)               v1(j) = yv(j, 1) + y_d_v(j, 1)
446               tair1(j) = yt(j, 1) + y_d_t(j, 1)               tair1(j) = yt(j, 1) + y_d_t(j, 1)
447               qair1(j) = yq(j, 1) + y_d_q(j, 1)               qair1(j) = yq(j, 1) + y_d_q(j, 1)
448               zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &               zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
449                    1)))*(ypaprs(j, 1)-ypplay(j, 1))                    1))) * (ypaprs(j, 1)-ypplay(j, 1))
450               tairsol(j) = yts(j) + y_d_ts(j)               tairsol(j) = yts(j) + y_d_ts(j)
451               rugo1(j) = yrugos(j)               rugo1(j) = yrugos(j)
452               IF (nsrf == is_oce) THEN               IF (nsrf == is_oce) THEN
453                  rugo1(j) = rugos(i, nsrf)                  rugo1(j) = frugs(i, nsrf)
454               END IF               END IF
455               psfce(j) = ypaprs(j, 1)               psfce(j) = ypaprs(j, 1)
456               patm(j) = ypplay(j, 1)               patm(j) = ypplay(j, 1)
# Line 681  contains Line 458  contains
458               qairsol(j) = yqsurf(j)               qairsol(j) = yqsurf(j)
459            END DO            END DO
460    
461            CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &            CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), qair1, &
462                 zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &                 zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, &
463                 yt10m, yq10m, yu10m, yustar)                 yq10m, wind10m(:knon), ustar(:knon))
464    
465            DO j = 1, knon            DO j = 1, knon
466               i = ni(j)               i = ni(j)
467               t2m(i, nsrf) = yt2m(j)               t2m(i, nsrf) = yt2m(j)
468               q2m(i, nsrf) = yq2m(j)               q2m(i, nsrf) = yq2m(j)
469    
470               ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman               u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
471               u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)                    / sqrt(u1(j)**2 + v1(j)**2)
472               v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)               v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
473                      / sqrt(u1(j)**2 + v1(j)**2)
474            END DO            END DO
475    
476            CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
477                 y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &                 y_flux_q(:knon), yu, yv, yt, yq, ypblh(:knon), ycapcl, &
478                 ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)                 yoliqcl, ycteicl, ypblt, ytherm, ylcl)
479    
480            DO j = 1, knon            DO j = 1, knon
481               i = ni(j)               i = ni(j)
# Line 709  contains Line 486  contains
486               cteicl(i, nsrf) = ycteicl(j)               cteicl(i, nsrf) = ycteicl(j)
487               pblt(i, nsrf) = ypblt(j)               pblt(i, nsrf) = ypblt(j)
488               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)  
489            END DO            END DO
490    
491            DO j = 1, knon            DO j = 1, knon
# Line 720  contains Line 494  contains
494                  q2(i, k, nsrf) = yq2(j, k)                  q2(i, k, nsrf) = yq2(j, k)
495               END DO               END DO
496            END DO            END DO
497            !IM "slab" ocean         else
498            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  
499         end IF if_knon         end IF if_knon
500      END DO loop_surface      END DO loop_surface
501    
502      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
503        frugs(:, is_oce) = rugmer
504        pctsrf(:, is_oce) = pctsrf_new_oce
505        pctsrf(:, is_sic) = pctsrf_new_sic
506    
507      rugos(:, is_oce) = rugmer      firstcal = .false.
     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.72  
changed lines
  Added in v.276
  ViewVC Help
Powered by ViewVC 1.1.21