/[lmdze]/trunk/Sources/phylmd/clmain.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/clmain.f

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

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