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

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

  ViewVC Help
Powered by ViewVC 1.1.21