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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21