/[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/libf/phylmd/clmain.f90 revision 38 by guez, Thu Jan 6 17:52:19 2011 UTC trunk/Sources/phylmd/clmain.f revision 248 by guez, Fri Jan 5 16:40:13 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        ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
20      ! Tout ce qui a trait aux traceurs est dans phytrac maintenant.      ! de la couche limite pour les traceurs se fait avec "cltrac" et
21      ! Pour l'instant le calcul de la couche limite pour les traceurs      ! ne tient pas compte de la diff\'erentiation des sous-fractions
22      ! se fait avec cltrac et ne tient pas compte de la différentiation      ! de sol.
23      ! des sous-fractions de sol.  
24        use clcdrag_m, only: clcdrag
25      ! Pour pouvoir extraire les coefficients d'échanges et le vent      use clqh_m, only: clqh
26      ! dans la première couche, trois champs supplémentaires ont été créés :      use clvent_m, only: clvent
27      ! zcoefh, zu1 et zv1. Pour l'instant nous avons moyenné les valeurs      use coefkz_m, only: coefkz
28      ! de ces trois champs sur les 4 sous-surfaces du modèle. Dans l'avenir      use coefkzmin_m, only: coefkzmin
29      ! si les informations des sous-surfaces doivent être prises en compte      use coefkz2_m, only: coefkz2
30      ! il faudra sortir ces mêmes champs en leur ajoutant une dimension,      USE conf_gcm_m, ONLY: lmt_pas
31      ! c'est a dire nbsrf (nombre de sous-surfaces).      USE conf_phys_m, ONLY: iflag_pbl
32        USE dimphy, ONLY: klev, klon, zmasq
33      ! Auteur Z.X. Li (LMD/CNRS) date: 1993/08/18      USE dimsoil, ONLY: nsoilmx
34      ! Objet : interface de "couche limite" (diffusion verticale)      use hbtm_m, only: hbtm
35        USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
36      ! Arguments:      USE interfoce_lim_m, ONLY: interfoce_lim
37      ! dtime----input-R- interval du temps (secondes)      use stdlevvar_m, only: stdlevvar
38      ! itap-----input-I- numero du pas de temps      USE suphec_m, ONLY: rd, rg, rkappa
39      ! date0----input-R- jour initial      use time_phylmdz, only: itap
40      ! t--------input-R- temperature (K)      use ustarhb_m, only: ustarhb
41      ! q--------input-R- vapeur d'eau (kg/kg)      use yamada4_m, only: yamada4
42      ! u--------input-R- vitesse u  
43      ! v--------input-R- vitesse v      REAL, INTENT(IN):: dtime ! interval du temps (secondes)
44      ! ts-------input-R- temperature du sol (en Kelvin)  
45      ! paprs----input-R- pression a intercouche (Pa)      REAL, INTENT(inout):: pctsrf(klon, nbsrf)
46      ! pplay----input-R- pression au milieu de couche (Pa)      ! tableau des pourcentages de surface de chaque maille
47      ! radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2  
48      ! rlat-----input-R- latitude en degree      REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
49      ! rugos----input-R- longeur de rugosite (en m)      REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg / kg)
50      ! cufi-----input-R- resolution des mailles en x (m)      REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
51      ! cvfi-----input-R- resolution des mailles en y (m)      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    
83        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      ! d_u------output-R- le changement pour "u"  
87      ! d_v------output-R- le changement pour "v"      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
88      ! d_ts-----output-R- le changement pour "ts"      ! changement pour "u" et "v"
89      ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)  
90      !                    (orientation positive vers le bas)      REAL, intent(out):: d_ts(:, :) ! (klon, nbsrf) variation of ftsol
91      ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)  
92      ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal      REAL, intent(out):: flux_t(klon, nbsrf)
93      ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal      ! 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
     ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')  
     ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')  
109    
110      ! tslab-in/output-R temperature du slab ocean (en Kelvin)      REAL, intent(out):: coefh(:, 2:) ! (klon, 2:klev)
111      ! uniqmnt pour slab      ! Pour pouvoir extraire les coefficients d'\'echange, le champ
112        ! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
113      ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')      ! ce champ sur les quatre sous-surfaces du mod\`ele.
114      !cc  
115      ! ffonte----Flux thermique utilise pour fondre la neige      REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)
116      ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la  
117      !           hauteur de neige, en kg/m2/s      REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)
118      ! on rajoute en output yu1 et yv1 qui sont les vents dans      ! composantes du vent \`a 10m sans spirale d'Ekman
119      ! la premiere couche  
120      ! ces 4 variables sont maintenant traites dans phytrac      ! Ionela Musat. Cf. Anne Mathieu : planetary boundary layer, hbtm.
121      ! itr--------input-I- nombre de traceurs      ! Comme les autres diagnostics on cumule dans physiq ce qui permet
122      ! tr---------input-R- q. de traceurs      ! de sortir les grandeurs par sous-surface.
123      ! flux_surf--input-R- flux de traceurs a la surface      REAL pblh(klon, nbsrf) ! height of planetary boundary layer
124      ! d_tr-------output-R tendance de traceurs      REAL capcl(klon, nbsrf)
125      !IM cf. AM : PBL      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      ! trmb1-------deep_cape
131        REAL trmb2(klon, nbsrf)
132      ! trmb2--------inhibition      ! trmb2--------inhibition
133        REAL trmb3(klon, nbsrf)
134      ! trmb3-------Point Omega      ! trmb3-------Point Omega
135      ! Cape(klon)-------Cape du thermique      REAL plcl(klon, nbsrf)
     ! EauLiq(klon)-------Eau liqu integr du thermique  
     ! ctei(klon)-------Critere d'instab d'entrainmt des nuages de CL  
     ! lcl------- Niveau de condensation  
     ! pblh------- HCL  
     ! pblT------- T au nveau HCL  
   
     USE histcom, ONLY : histbeg_totreg, histdef, histend, histsync  
     use histwrite_m, only: histwrite  
     use calendar, ONLY : ymds2ju  
     USE dimens_m, ONLY : iim, jjm  
     USE indicesol, ONLY : epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf  
     USE dimphy, ONLY : klev, klon, zmasq  
     USE dimsoil, ONLY : nsoilmx  
     USE temps, ONLY : annee_ref, itau_phy  
     USE dynetat0_m, ONLY : day_ini  
     USE iniprint, ONLY : prt_level  
     USE suphec_m, ONLY : rd, rg, rkappa  
     USE conf_phys_m, ONLY : iflag_pbl  
     USE gath_cpl, ONLY : gath2cpl  
     use hbtm_m, only: hbtm  
   
     REAL, INTENT (IN) :: dtime  
     REAL date0  
     INTEGER, INTENT (IN) :: itap  
     REAL t(klon, klev), q(klon, klev)  
     REAL u(klon, klev), v(klon, klev)  
     REAL, INTENT (IN) :: paprs(klon, klev+1)  
     REAL, INTENT (IN) :: pplay(klon, klev)  
     REAL, INTENT (IN) :: rlon(klon), rlat(klon)  
     REAL cufi(klon), cvfi(klon)  
     REAL d_t(klon, klev), d_q(klon, klev)  
     REAL d_u(klon, klev), d_v(klon, klev)  
     REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)  
     REAL dflux_t(klon), dflux_q(klon)  
     !IM "slab" ocean  
     REAL flux_o(klon), flux_g(klon)  
     REAL y_flux_o(klon), y_flux_g(klon)  
     REAL tslab(klon), ytslab(klon)  
     REAL seaice(klon), y_seaice(klon)  
     REAL y_fqcalving(klon), y_ffonte(klon)  
136      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
137      REAL run_off_lic_0(klon), y_run_off_lic_0(klon)      ! ffonte----Flux thermique utilise pour fondre la neige
138        ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
139      REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)      !           hauteur de neige, en kg / m2 / s
140      REAL rugmer(klon), agesno(klon, nbsrf)      REAL run_off_lic_0(klon)
     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)  
     REAL d_ts(klon, nbsrf)  
     REAL snow(klon, nbsrf)  
     REAL qsurf(klon, nbsrf)  
     REAL evap(klon, nbsrf)  
     REAL albe(klon, nbsrf)  
     REAL alblw(klon, nbsrf)  
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)  
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)
     REAL zcoefh(klon, klev)  
     REAL zu1(klon)  
     REAL zv1(klon)  
   
     !$$$ PB ajout pour soil  
     LOGICAL, INTENT (IN) :: soil_model  
     !IM ajout seuils cdrm, cdrh  
     REAL cdmmax, cdhmax  
149    
150      REAL ksta, ksta_ter      REAL y_fqcalving(klon), y_ffonte(klon)
151      LOGICAL ok_kzmin      real y_run_off_lic_0(klon)
152        REAL rugmer(klon)
     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 clqh, clvent, coefkz, 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      REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)      REAL yrain_f(klon) ! liquid water mass flux (kg / m2 / s), positive down
160      REAL yrain_f(klon), ysnow_f(klon)      REAL ysnow_f(klon) ! solid water mass flux (kg / m2 / s), positive down
     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)  
   
     !IM 081204 hcl_Anne ? BEG  
     REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)  
     REAL ykmm(klon, klev+1), ykmn(klon, klev+1)  
     REAL ykmq(klon, klev+1)  
     REAL yq2(klon, klev+1), q2(klon, klev+1, nbsrf)  
     REAL q2diag(klon, klev+1)  
     !IM 081204 hcl_Anne ? END  
   
     REAL u1lay(klon), v1lay(klon)  
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
     ! Introduction d'une variable "pourcentage potentiel" pour tenir compte  
     ! des eventuelles apparitions et/ou disparitions de la glace de mer  
     REAL pctsrf_pot(klon, nbsrf)  
180    
181      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.      REAL pctsrf_pot(klon, nbsrf)
182        ! "pourcentage potentiel" pour tenir compte des \'eventuelles
183      ! maf pour sorties IOISPL en cas de debugagage      ! apparitions ou disparitions de la glace de mer
184    
185      CHARACTER (80) cldebug      REAL yt2m(klon), yq2m(klon), wind10m(klon)
186      SAVE cldebug      REAL ustar(klon)
     CHARACTER (8) cl_surf(nbsrf)  
     SAVE cl_surf  
     INTEGER nhoridbg, nidbg  
     SAVE nhoridbg, nidbg  
     INTEGER ndexbg(iim*(jjm+1))  
     REAL zx_lon(iim, jjm+1), zx_lat(iim, jjm+1), zjulian  
     REAL tabindx(klon)  
     REAL debugtab(iim, jjm+1)  
     LOGICAL first_appel  
     SAVE first_appel  
     DATA first_appel/ .TRUE./  
     LOGICAL :: debugindex = .FALSE.  
     INTEGER idayref  
     REAL t2m(klon, nbsrf), q2m(klon, nbsrf)  
     REAL u10m(klon, nbsrf), v10m(klon, nbsrf)  
   
     REAL yt2m(klon), yq2m(klon), yu10m(klon)  
     REAL yustar(klon)  
     ! -- LOOP  
     REAL yu10mx(klon)  
     REAL yu10my(klon)  
     REAL ywindsp(klon)  
     ! -- LOOP  
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)  
     REAL plcl(klon, nbsrf)  
     REAL capcl(klon, nbsrf)  
     REAL oliqcl(klon, nbsrf)  
     REAL cteicl(klon, nbsrf)  
     REAL pblt(klon, nbsrf)  
     REAL therm(klon, nbsrf)  
     REAL trmb1(klon, nbsrf)  
     REAL trmb2(klon, nbsrf)  
     REAL trmb3(klon, nbsrf)  
189      REAL ypblh(klon)      REAL ypblh(klon)
190      REAL ylcl(klon)      REAL ylcl(klon)
191      REAL ycapcl(klon)      REAL ycapcl(klon)
# Line 278  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    
     ! initialisation Anne  
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.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  
   
     ! initialisation:  
216    
217      DO i = 1, klon      ! Initialization:
218         rugmer(i) = 0.0      rugmer = 0.
219         cdragh(i) = 0.0      cdragh = 0.
220         cdragm(i) = 0.0      cdragm = 0.
221         dflux_t(i) = 0.0      dflux_t = 0.
222         dflux_q(i) = 0.0      dflux_q = 0.
223         zu1(i) = 0.0      ypct = 0.
224         zv1(i) = 0.0      yqsurf = 0.
225      END DO      yrain_f = 0.
226      ypct = 0.0      ysnow_f = 0.
227      yts = 0.0      yrugos = 0.
228      ysnow = 0.0      ypaprs = 0.
229      yqsurf = 0.0      ypplay = 0.
230      yalb = 0.0      ydelp = 0.
231      yalblw = 0.0      yu = 0.
232      yrain_f = 0.0      yv = 0.
233      ysnow_f = 0.0      yt = 0.
234      yfder = 0.0      yq = 0.
235      ytaux = 0.0      y_dflux_t = 0.
236      ytauy = 0.0      y_dflux_q = 0.
     ysolsw = 0.0  
     ysollw = 0.0  
     ysollwdown = 0.0  
     yrugos = 0.0  
     yu1 = 0.0  
     yv1 = 0.0  
     yrads = 0.0  
     ypaprs = 0.0  
     ypplay = 0.0  
     ydelp = 0.0  
     yu = 0.0  
     yv = 0.0  
     yt = 0.0  
     yq = 0.0  
     pctsrf_new = 0.0  
     y_flux_u = 0.0  
     y_flux_v = 0.0  
     !$$ PB  
     y_dflux_t = 0.0  
     y_dflux_q = 0.0  
     ytsoil = 999999.  
237      yrugoro = 0.      yrugoro = 0.
238      ! -- LOOP      d_ts = 0.
     yu10mx = 0.0  
     yu10my = 0.0  
     ywindsp = 0.0  
     ! -- LOOP  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           d_ts(i, nsrf) = 0.0  
        END DO  
     END DO  
     !§§§ PB  
     yfluxlat = 0.  
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      DO k = 1, klev      fluxlat = 0.
244         DO i = 1, klon      d_t = 0.
245            d_t(i, k) = 0.0      d_q = 0.
246            d_q(i, k) = 0.0      d_u = 0.
247            d_u(i, k) = 0.0      d_v = 0.
248            d_v(i, k) = 0.0      coefh = 0.
249            zcoefh(i, k) = 0.0  
250         END DO      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
251      END DO      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
252        ! (\`a affiner)
     ! Boucler sur toutes les sous-fractions du sol:  
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      DO nsrf = 1, nbsrf      ! Tester si c'est le moment de lire le fichier:
260         ! chercher les indices:      if (mod(itap - 1, lmt_pas) == 0) then
261           CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
262        endif
263    
264        ! Boucler sur toutes les sous-fractions du sol:
265    
266        loop_surface: DO nsrf = 1, nbsrf
267           ! Chercher les indices :
268         ni = 0         ni = 0
269         knon = 0         knon = 0
270         DO i = 1, klon         DO i = 1, klon
271            ! pour determiner le domaine a 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 435  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  
           DO j = 1, knon  
              i = ni(j)  
              yqsol(j) = qsol(i)  
           END DO  
        ELSE  
           yqsol = 0.  
        END IF  
        !$$$ PB ajour pour soil  
        DO k = 1, nsoilmx  
           DO j = 1, knon  
              i = ni(j)  
              ytsoil(j, k) = ftsoil(i, k, nsrf)  
           END DO  
        END DO  
        DO k = 1, klev  
280            DO j = 1, knon            DO j = 1, knon
281               i = ni(j)               i = ni(j)
282               ypaprs(j, k) = paprs(i, k)               ypct(j) = pctsrf(i, nsrf)
283               ypplay(j, k) = pplay(i, k)               yts(j) = ftsol(i, nsrf)
284               ydelp(j, k) = delp(i, k)               snow(j) = fsnow(i, nsrf)
285               yu(j, k) = u(i, k)               yqsurf(j) = qsurf(i, nsrf)
286               yv(j, k) = v(i, k)               yalb(j) = falbe(i, nsrf)
287               yt(j, k) = t(i, k)               yrain_f(j) = rain_fall(i)
288               yq(j, k) = q(i, k)               ysnow_f(j) = snow_f(i)
289                 yagesno(j) = agesno(i, nsrf)
290                 yrugos(j) = frugs(i, nsrf)
291                 yrugoro(j) = rugoro(i)
292                 yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)
293                 ypaprs(j, klev + 1) = paprs(i, klev + 1)
294                 y_run_off_lic_0(j) = run_off_lic_0(i)
295            END DO            END DO
        END DO  
296    
297         ! calculer Cdrag et les coefficients d'echange            ! For continent, copy soil water content
298         CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,&            IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
             yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)  
        !IM 081204 BEG  
        !CR test  
        IF (iflag_pbl==1) THEN  
           !IM 081204 END  
           CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)  
           DO k = 1, klev  
              DO i = 1, knon  
                 ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))  
                 ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))  
              END DO  
           END DO  
        END IF  
299    
300         !IM cf JLD : on seuille ycoefm et ycoefh            ytsoil(:knon, :) = ftsoil(ni(:knon), :, nsrf)
        IF (nsrf==is_oce) THEN  
           DO j = 1, knon  
              !           ycoefm(j, 1)=min(ycoefm(j, 1), 1.1E-3)  
              ycoefm(j, 1) = min(ycoefm(j, 1), cdmmax)  
              !           ycoefh(j, 1)=min(ycoefh(j, 1), 1.1E-3)  
              ycoefh(j, 1) = min(ycoefh(j, 1), cdhmax)  
           END DO  
        END IF  
   
        !IM: 261103  
        IF (ok_kzmin) THEN  
           !IM cf FH: 201103 BEG  
           !   Calcul d'une diffusion minimale pour les conditions tres stables.  
           CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, ycoefm, &  
                ycoefm0, ycoefh0)  
301    
           IF (1==1) THEN  
              DO k = 1, klev  
                 DO i = 1, knon  
                    ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))  
                    ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))  
                 END DO  
              END DO  
           END IF  
           !IM cf FH: 201103 END  
           !IM: 261103  
        END IF !ok_kzmin  
   
        IF (iflag_pbl>=3) THEN  
           ! MELLOR ET YAMADA adapté à Mars, Richard Fournier et Frédéric Hourdin  
           yzlay(1:knon, 1) = rd*yt(1:knon, 1)/(0.5*(ypaprs(1:knon, &  
                1)+ypplay(1:knon, 1)))*(ypaprs(1:knon, 1)-ypplay(1:knon, 1))/rg  
           DO k = 2, klev  
              yzlay(1:knon, k) = yzlay(1:knon, k-1) &  
                   + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &  
                   / ypaprs(1:knon, k) &  
                   * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg  
           END DO  
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            !   Bug introduit volontairement pour converger avec les resultats            ! Calculer les géopotentiels de chaque couche:
316            !  du papier sur les thermiques.  
317            IF (1==1) THEN            zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
318               y_cd_m(1:knon) = ycoefm(1:knon, 1)                 + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
              y_cd_h(1:knon) = ycoefh(1:knon, 1)  
           ELSE  
              y_cd_h(1:knon) = ycoefm(1:knon, 1)  
              y_cd_m(1:knon) = ycoefh(1:knon, 1)  
           END IF  
           CALL ustarhb(knon, yu, yv, y_cd_m, yustar)  
319    
320            IF (prt_level>9) THEN            DO k = 2, klev
321               PRINT *, 'USTAR = ', yustar               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              CALL coefkz(nsrf, ypaprs(:knon, :), ypplay(:knon, :), ksta, &
331                   ksta_ter, yts(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
332                   yq(:knon, :), zgeop(:knon, :), ycoefm(:knon, :), &
333                   ycoefh(:knon, :))
334    
335              IF (iflag_pbl == 1) THEN
336                 CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0(:knon, :), &
337                      ycoefh0(:knon, :))
338                 ycoefm(:knon, :) = max(ycoefm(:knon, :), ycoefm0(:knon, :))
339                 ycoefh(:knon, :) = max(ycoefh(:knon, :), ycoefh0(:knon, :))
340                 ycdragm(:knon) = max(ycdragm(:knon), 0.)
341                 ycdragh(:knon) = max(ycdragh(:knon), 0.)
342            END IF            END IF
343    
344            !   iflag_pbl peut etre utilise comme longuer de melange            ! on met un seuil pour ycdragm et ycdragh
345              IF (nsrf == is_oce) THEN
346                 ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
347                 ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
348              END IF
349    
350            IF (iflag_pbl>=11) THEN            IF (ok_kzmin) THEN
351               CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &               ! Calcul d'une diffusion minimale pour les conditions tres stables
352                    yu, yv, yteta, y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, &               CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
353                    iflag_pbl)                    ycdragm(:knon), ycoefh0(:knon, :))
354            ELSE               ycoefm0(:knon, :) = ycoefh0(:knon, :)
355               CALL yamada4(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, &               ycoefm(:knon, :) = max(ycoefm(:knon, :), ycoefm0(:knon, :))
356                    yv, yteta, y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)               ycoefh(:knon, :) = max(ycoefh(:knon, :), ycoefh0(:knon, :))
357            END IF            END IF
358    
359            ycoefm(1:knon, 1) = y_cd_m(1:knon)            IF (iflag_pbl >= 6) THEN
360            ycoefh(1:knon, 1) = y_cd_h(1:knon)               ! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et
361            ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)               ! Fr\'ed\'eric Hourdin
362            ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)               yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
363         END IF                    + ypplay(:knon, 1))) &
364                      * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
365         ! calculer la diffusion des vitesses "u" et "v"  
366         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &               DO k = 2, klev
367              ydelp, y_d_u, y_flux_u)                  yzlay(:knon, k) = yzlay(:knon, k-1) &
368         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &                       + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
369              ydelp, y_d_v, y_flux_v)                       / ypaprs(1:knon, k) &
370                         * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
371         ! pour le couplage               END DO
        ytaux = y_flux_u(:, 1)  
        ytauy = y_flux_v(:, 1)  
   
        ! FH modif sur le cdrag temperature  
        !$$$PB : déplace dans clcdrag  
        !$$$      do i=1, knon  
        !$$$         ycoefh(i, 1)=ycoefm(i, 1)*0.8  
        !$$$      enddo  
   
        ! 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  
372    
373         DO k = 1, klev               DO k = 1, klev
374            DO j = 1, knon                  yteta(1:knon, k) = yt(1:knon, k) * (ypaprs(1:knon, 1) &
375               i = ni(j)                       / ypplay(1:knon, k))**rkappa * (1. + 0.61 * yq(1:knon, k))
376               ycoefh(j, k) = ycoefh(j, k)*ypct(j)               END DO
              ycoefm(j, k) = ycoefm(j, k)*ypct(j)  
              y_d_t(j, k) = y_d_t(j, k)*ypct(j)  
              y_d_q(j, k) = y_d_q(j, k)*ypct(j)  
              !§§§ PB  
              flux_t(i, k, nsrf) = y_flux_t(j, k)  
              flux_q(i, k, nsrf) = y_flux_q(j, k)  
              flux_u(i, k, nsrf) = y_flux_u(j, k)  
              flux_v(i, k, nsrf) = y_flux_v(j, k)  
              !$$$ PB        y_flux_t(j, k) = y_flux_t(j, k) * ypct(j)  
              !$$$ PB        y_flux_q(j, k) = y_flux_q(j, k) * ypct(j)  
              y_d_u(j, k) = y_d_u(j, k)*ypct(j)  
              y_d_v(j, k) = y_d_v(j, k)*ypct(j)  
              !$$$ PB        y_flux_u(j, k) = y_flux_u(j, k) * ypct(j)  
              !$$$ PB        y_flux_v(j, k) = y_flux_v(j, k) * ypct(j)  
           END DO  
        END DO  
377    
378         evap(:, nsrf) = -flux_q(:, 1, nsrf)               zlev(:knon, 1) = 0.
379                 zlev(:knon, klev + 1) = 2. * yzlay(:knon, klev) &
380                      - yzlay(:knon, klev - 1)
381    
382         albe(:, nsrf) = 0.               DO k = 2, klev
383         alblw(:, nsrf) = 0.                  zlev(:knon, k) = 0.5 * (yzlay(:knon, k) + yzlay(:knon, k-1))
384         snow(:, nsrf) = 0.               END DO
385         qsurf(:, nsrf) = 0.  
386         rugos(:, nsrf) = 0.               DO k = 1, klev + 1
387         fluxlat(:, nsrf) = 0.                  DO j = 1, knon
388         DO j = 1, knon                     i = ni(j)
389            i = ni(j)                     yq2(j, k) = q2(i, k, nsrf)
390            d_ts(i, nsrf) = y_d_ts(j)                  END DO
391            albe(i, nsrf) = yalb(j)               END DO
392            alblw(i, nsrf) = yalblw(j)  
393            snow(i, nsrf) = ysnow(j)               ustar(:knon) = ustarhb(yu(:knon, 1), yv(:knon, 1), ycdragm(:knon))
394            qsurf(i, nsrf) = yqsurf(j)               CALL yamada4(dtime, rg, zlev(:knon, :), yzlay(:knon, :), &
395            rugos(i, nsrf) = yz0_new(j)                    yu(:knon, :), yv(:knon, :), yteta(:knon, :), yq2(:knon, :), &
396            fluxlat(i, nsrf) = yfluxlat(j)                    ycoefm(:knon, :), ycoefh(:knon, :), ustar(:knon))
397            !$$$ pb         rugmer(i) = yrugm(j)            END IF
398            IF (nsrf==is_oce) THEN  
399               rugmer(i) = yrugm(j)            CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
400               rugos(i, nsrf) = yrugm(j)                 ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
401                   ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
402                   y_flux_u(:knon))
403              CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
404                   ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
405                   ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
406                   y_flux_v(:knon))
407    
408              ! calculer la diffusion de "q" et de "h"
409              CALL clqh(dtime, julien, firstcal, nsrf, ni(:knon), &
410                   ytsoil(:knon, :), yqsol(:knon), mu0, yrugos, yrugoro, &
411                   yu(:knon, 1), yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), &
412                   yt, yq, yts(:knon), ypaprs, ypplay, ydelp, yrads(:knon), &
413                   yalb(:knon), snow(:knon), yqsurf, yrain_f, ysnow_f, &
414                   yfluxlat(:knon), pctsrf_new_sic, yagesno(:knon), y_d_t, y_d_q, &
415                   y_d_ts(:knon), yz0_new, y_flux_t(:knon), y_flux_q(:knon), &
416                   y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving, y_ffonte, &
417                   y_run_off_lic_0)
418    
419              ! calculer la longueur de rugosite sur ocean
420              yrugm = 0.
421              IF (nsrf == is_oce) THEN
422                 DO j = 1, knon
423                    yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
424                         / rg + 0.11 * 14E-6 &
425                         / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
426                    yrugm(j) = max(1.5E-05, yrugm(j))
427                 END DO
428            END IF            END IF
           !IM cf JLD ??  
           agesno(i, nsrf) = yagesno(j)  
           fqcalving(i, nsrf) = y_fqcalving(j)  
           ffonte(i, nsrf) = y_ffonte(j)  
           cdragh(i) = cdragh(i) + ycoefh(j, 1)  
           cdragm(i) = cdragm(i) + ycoefm(j, 1)  
           dflux_t(i) = dflux_t(i) + y_dflux_t(j)  
           dflux_q(i) = dflux_q(i) + y_dflux_q(j)  
           zu1(i) = zu1(i) + yu1(j)  
           zv1(i) = zv1(i) + yv1(j)  
        END DO  
        IF (nsrf==is_ter) THEN  
           DO j = 1, knon  
              i = ni(j)  
              qsol(i) = yqsol(j)  
           END DO  
        END IF  
        IF (nsrf==is_lic) THEN  
           DO j = 1, knon  
              i = ni(j)  
              run_off_lic_0(i) = y_run_off_lic_0(j)  
           END DO  
        END IF  
        !$$$ PB ajout pour soil  
        ftsoil(:, :, nsrf) = 0.  
        DO k = 1, nsoilmx  
429            DO j = 1, knon            DO j = 1, knon
430               i = ni(j)               y_dflux_t(j) = y_dflux_t(j) * ypct(j)
431               ftsoil(i, k, nsrf) = ytsoil(j, k)               y_dflux_q(j) = y_dflux_q(j) * ypct(j)
432            END DO            END DO
        END DO  
433    
        DO j = 1, knon  
           i = ni(j)  
434            DO k = 1, klev            DO k = 1, klev
435               d_t(i, k) = d_t(i, k) + y_d_t(j, k)               DO j = 1, knon
436               d_q(i, k) = d_q(i, k) + y_d_q(j, k)                  i = ni(j)
437               !$$$ PB        flux_t(i, k) = flux_t(i, k) + y_flux_t(j, k)                  y_d_t(j, k) = y_d_t(j, k) * ypct(j)
438               !$$$         flux_q(i, k) = flux_q(i, k) + y_flux_q(j, k)                  y_d_q(j, k) = y_d_q(j, k) * ypct(j)
439               d_u(i, k) = d_u(i, k) + y_d_u(j, k)                  y_d_u(j, k) = y_d_u(j, k) * ypct(j)
440               d_v(i, k) = d_v(i, k) + y_d_v(j, k)                  y_d_v(j, k) = y_d_v(j, k) * ypct(j)
441               !$$$  PB       flux_u(i, k) = flux_u(i, k) + y_flux_u(j, k)               END DO
              !$$$         flux_v(i, k) = flux_v(i, k) + y_flux_v(j, k)  
              zcoefh(i, k) = zcoefh(i, k) + ycoefh(j, k)  
442            END DO            END DO
        END DO  
443    
444         !cc diagnostic t, q a 2m et u, v a 10m            flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
445              flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
446         DO j = 1, knon            flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
447            i = ni(j)            flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
448            uzon(j) = yu(j, 1) + y_d_u(j, 1)  
449            vmer(j) = yv(j, 1) + y_d_v(j, 1)            evap(:, nsrf) = -flux_q(:, nsrf)
450            tair1(j) = yt(j, 1) + y_d_t(j, 1)  
451            qair1(j) = yq(j, 1) + y_d_q(j, 1)            falbe(:, nsrf) = 0.
452            zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &            fsnow(:, nsrf) = 0.
453                 1)))*(ypaprs(j, 1)-ypplay(j, 1))            qsurf(:, nsrf) = 0.
454            tairsol(j) = yts(j) + y_d_ts(j)            frugs(:, nsrf) = 0.
455            rugo1(j) = yrugos(j)            DO j = 1, knon
456            IF (nsrf==is_oce) THEN               i = ni(j)
457               rugo1(j) = rugos(i, nsrf)               d_ts(i, nsrf) = y_d_ts(j)
458                 falbe(i, nsrf) = yalb(j)
459                 fsnow(i, nsrf) = snow(j)
460                 qsurf(i, nsrf) = yqsurf(j)
461                 frugs(i, nsrf) = yz0_new(j)
462                 fluxlat(i, nsrf) = yfluxlat(j)
463                 IF (nsrf == is_oce) THEN
464                    rugmer(i) = yrugm(j)
465                    frugs(i, nsrf) = yrugm(j)
466                 END IF
467                 agesno(i, nsrf) = yagesno(j)
468                 fqcalving(i, nsrf) = y_fqcalving(j)
469                 ffonte(i, nsrf) = y_ffonte(j)
470                 cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
471                 cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
472                 dflux_t(i) = dflux_t(i) + y_dflux_t(j)
473                 dflux_q(i) = dflux_q(i) + y_dflux_q(j)
474              END DO
475              IF (nsrf == is_ter) THEN
476                 qsol(ni(:knon)) = yqsol(:knon)
477              else IF (nsrf == is_lic) THEN
478                 DO j = 1, knon
479                    i = ni(j)
480                    run_off_lic_0(i) = y_run_off_lic_0(j)
481                 END DO
482            END IF            END IF
           psfce(j) = ypaprs(j, 1)  
           patm(j) = ypplay(j, 1)  
483    
484            qairsol(j) = yqsurf(j)            ftsoil(:, :, nsrf) = 0.
485         END DO            ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)
486    
487         CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &            DO j = 1, knon
488              tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &               i = ni(j)
489              yu10m, yustar)               DO k = 1, klev
490         !IM 081204 END                  d_t(i, k) = d_t(i, k) + y_d_t(j, k)
491                    d_q(i, k) = d_q(i, k) + y_d_q(j, k)
492         DO j = 1, knon                  d_u(i, k) = d_u(i, k) + y_d_u(j, k)
493            i = ni(j)                  d_v(i, k) = d_v(i, k) + y_d_v(j, k)
494            t2m(i, nsrf) = yt2m(j)               END DO
495            q2m(i, nsrf) = yq2m(j)            END DO
   
           ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman  
           u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)  
           v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)  
496    
497         END DO            forall (k = 2:klev) coefh(ni(:knon), k) &
498                   = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypct(:knon)
499    
500         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  
501    
502         DO j = 1, knon            DO j = 1, knon
           DO k = 1, klev + 1  
503               i = ni(j)               i = ni(j)
504               q2(i, k, nsrf) = yq2(j, k)               u1(j) = yu(j, 1) + y_d_u(j, 1)
505                 v1(j) = yv(j, 1) + y_d_v(j, 1)
506                 tair1(j) = yt(j, 1) + y_d_t(j, 1)
507                 qair1(j) = yq(j, 1) + y_d_q(j, 1)
508                 zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
509                      1))) * (ypaprs(j, 1)-ypplay(j, 1))
510                 tairsol(j) = yts(j) + y_d_ts(j)
511                 rugo1(j) = yrugos(j)
512                 IF (nsrf == is_oce) THEN
513                    rugo1(j) = frugs(i, nsrf)
514                 END IF
515                 psfce(j) = ypaprs(j, 1)
516                 patm(j) = ypplay(j, 1)
517    
518                 qairsol(j) = yqsurf(j)
519            END DO            END DO
520         END DO  
521         !IM "slab" ocean            CALL stdlevvar(klon, knon, nsrf, u1(:knon), v1(:knon), tair1(:knon), &
522         IF (nsrf==is_oce) THEN                 qair1, zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, &
523                   yq2m, yt10m, yq10m, wind10m(:knon), ustar(:knon))
524    
525            DO j = 1, knon            DO j = 1, knon
              ! on projette sur la grille globale  
526               i = ni(j)               i = ni(j)
527               IF (pctsrf_new(i, is_oce)>epsfra) THEN               t2m(i, nsrf) = yt2m(j)
528                  flux_o(i) = y_flux_o(j)               q2m(i, nsrf) = yq2m(j)
529               ELSE  
530                  flux_o(i) = 0.               u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
531               END IF                    / sqrt(u1(j)**2 + v1(j)**2)
532                 v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
533                      / sqrt(u1(j)**2 + v1(j)**2)
534            END DO            END DO
        END IF  
535    
536         IF (nsrf==is_sic) THEN            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
537                   y_flux_q(:knon), yu, yv, yt, yq, ypblh(:knon), ycapcl, &
538                   yoliqcl, ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
539    
540            DO j = 1, knon            DO j = 1, knon
541               i = ni(j)               i = ni(j)
542               ! On pondère lorsque l'on fait le bilan au sol :               pblh(i, nsrf) = ypblh(j)
543               ! flux_g(i) = y_flux_g(j)*ypct(j)               plcl(i, nsrf) = ylcl(j)
544               IF (pctsrf_new(i, is_sic)>epsfra) THEN               capcl(i, nsrf) = ycapcl(j)
545                  flux_g(i) = y_flux_g(j)               oliqcl(i, nsrf) = yoliqcl(j)
546               ELSE               cteicl(i, nsrf) = ycteicl(j)
547                  flux_g(i) = 0.               pblt(i, nsrf) = ypblt(j)
548               END IF               therm(i, nsrf) = ytherm(j)
549                 trmb1(i, nsrf) = ytrmb1(j)
550                 trmb2(i, nsrf) = ytrmb2(j)
551                 trmb3(i, nsrf) = ytrmb3(j)
552            END DO            END DO
553    
554         END IF            DO j = 1, knon
555         !nsrf.EQ.is_sic                                                           DO k = 1, klev + 1
556         IF (ocean=='slab  ') THEN                  i = ni(j)
557            IF (nsrf==is_oce) THEN                  q2(i, k, nsrf) = yq2(j, k)
558               tslab(1:klon) = ytslab(1:klon)               END DO
559               seaice(1:klon) = y_seaice(1:klon)            END DO
560               !nsrf                                                               else
561            END IF            fsnow(:, nsrf) = 0.
562            !OCEAN                                                               end IF if_knon
563         END IF      END DO loop_surface
     END DO  
564    
565      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
566      ! A rajouter: conservation de l'albedo      frugs(:, is_oce) = rugmer
567        pctsrf(:, is_oce) = pctsrf_new_oce
568        pctsrf(:, is_sic) = pctsrf_new_sic
569    
570      rugos(:, is_oce) = rugmer      firstcal = .false.
     pctsrf = pctsrf_new  
571    
572    END SUBROUTINE clmain    END SUBROUTINE clmain
573    

Legend:
Removed from v.38  
changed lines
  Added in v.248

  ViewVC Help
Powered by ViewVC 1.1.21