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

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

  ViewVC Help
Powered by ViewVC 1.1.21