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

Legend:
Removed from v.40  
changed lines
  Added in v.206

  ViewVC Help
Powered by ViewVC 1.1.21