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

Diff of /trunk/phylmd/pbl_surface.f

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

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

Legend:
Removed from v.57  
changed lines
  Added in v.249

  ViewVC Help
Powered by ViewVC 1.1.21