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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21