/[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

revision 40 by guez, Tue Feb 22 13:49:36 2011 UTC revision 62 by guez, Thu Jul 26 14:37:37 2012 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, date0, pctsrf, pctsrf_new, t, q, u, v, &
8         jour, rmu0, co2_ppm, ok_veget, ocean, npas, nexca, ts,&         jour, rmu0, co2_ppm, ok_veget, ocean, npas, nexca, ts, &
9         soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil,&         soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, &
10         qsol, paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat,&         qsol, paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat, &
11         rain_f, snow_f, solsw, sollw, sollwdown, fder, rlon, rlat, cufi,&         rain_fall, snow_f, solsw, sollw, sollwdown, fder, rlon, rlat, cufi, &
12         cvfi, rugos, debut, lafin, agesno, rugoro, d_t, d_q, d_u, d_v,&         cvfi, rugos, debut, lafin, agesno, rugoro, d_t, d_q, d_u, d_v, &
13         d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2,&         d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, &
14         dflux_t, dflux_q, zcoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh,&         dflux_t, dflux_q, zcoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh, &
15         capcl, oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl,&         capcl, oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl, &
16         fqcalving, ffonte, run_off_lic_0, flux_o, flux_g, tslab, seaice)         fqcalving, ffonte, run_off_lic_0, flux_o, flux_g, tslab, seaice)
17    
18      ! From phylmd/clmain.F, version 1.6 2005/11/16 14:47:19      ! From phylmd/clmain.F, version 1.6 2005/11/16 14:47:19
19      ! Author: Z.X. Li (LMD/CNRS), date: 1993/08/18      ! Author: Z. X. Li (LMD/CNRS), date: 1993/08/18
20      ! Objet : interface de "couche limite" (diffusion verticale)      ! Objet : interface de couche limite (diffusion verticale)
21    
22      ! Tout ce qui a trait aux traceurs est dans "phytrac" maintenant.      ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
23      ! Pour l'instant le calcul de la couche limite pour les traceurs      ! de la couche limite pour les traceurs se fait avec "cltrac" et
24      ! se fait avec "cltrac" et ne tient pas compte de la différentiation      ! ne tient pas compte de la différentiation des sous-fractions de
25      ! des sous-fractions de sol.      ! sol.
26    
27      ! Pour pouvoir extraire les coefficients d'échanges et le vent      ! Pour pouvoir extraire les coefficients d'échanges et le vent
28      ! dans la première couche, trois champs supplémentaires ont été      ! dans la première couche, trois champs ont été créés : "zcoefh",
29      ! créés : "zcoefh", "zu1" et "zv1". Pour l'instant nous avons      ! "zu1" et "zv1". Nous avons moyenné les valeurs de ces trois
30      ! moyenné les valeurs de ces trois champs sur les 4 sous-surfaces      ! champs sur les quatre sous-surfaces du modèle.
31      ! du modèle. Dans l'avenir, si les informations des sous-surfaces  
32      ! doivent être prises en compte, il faudra sortir ces mêmes champs      use calendar, ONLY: ymds2ju
33      ! en leur ajoutant une dimension, c'est-à-dire "nbsrf" (nombre de      use clqh_m, only: clqh
34      ! sous-surfaces).      use clvent_m, only: clvent
35        use coefkz_m, only: coefkz
36        use coefkzmin_m, only: coefkzmin
37        USE conf_gcm_m, ONLY: prt_level
38        USE conf_phys_m, ONLY: iflag_pbl
39        USE dimens_m, ONLY: iim, jjm
40        USE dimphy, ONLY: klev, klon, zmasq
41        USE dimsoil, ONLY: nsoilmx
42        USE dynetat0_m, ONLY: day_ini
43        USE gath_cpl, ONLY: gath2cpl
44        use hbtm_m, only: hbtm
45        USE histbeg_totreg_m, ONLY: histbeg_totreg
46        USE histdef_m, ONLY: histdef
47        USE histend_m, ONLY: histend
48        USE histsync_m, ONLY: histsync
49        use histwrite_m, only: histwrite
50        USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
51        USE suphec_m, ONLY: rd, rg, rkappa
52        USE temps, ONLY: annee_ref, itau_phy
53        use ustarhb_m, only: ustarhb
54        use vdif_kcay_m, only: vdif_kcay
55        use yamada4_m, only: yamada4
56    
57      ! Arguments:      ! Arguments:
58      ! dtime----input-R- interval du temps (secondes)  
59      ! itap-----input-I- numero du pas de temps      REAL, INTENT(IN):: dtime ! interval du temps (secondes)
60      ! date0----input-R- jour initial      INTEGER, INTENT(IN):: itap ! numero du pas de temps
61      ! t--------input-R- temperature (K)      REAL, INTENT(IN):: date0 ! jour initial
62      ! q--------input-R- vapeur d'eau (kg/kg)      REAL, INTENT(inout):: pctsrf(klon, nbsrf)
63      ! u--------input-R- vitesse u  
64      ! v--------input-R- vitesse v      ! la nouvelle repartition des surfaces sortie de l'interface
65      ! ts-------input-R- temperature du sol (en Kelvin)      REAL, INTENT(out):: pctsrf_new(klon, nbsrf)
66      ! paprs----input-R- pression a intercouche (Pa)  
67      ! pplay----input-R- pression au milieu de couche (Pa)      REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
68      ! radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2      REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg/kg)
69      ! rlat-----input-R- latitude en degree      REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
70      ! rugos----input-R- longeur de rugosite (en m)      INTEGER, INTENT(IN):: jour ! jour de l'annee en cours
71        REAL, intent(in):: rmu0(klon) ! cosinus de l'angle solaire zenithal    
72        REAL, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)
73        REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
74        REAL, INTENT(IN):: rlon(klon)
75        REAL, INTENT(IN):: rlat(klon) ! latitude en degrés
76        REAL cufi(klon), cvfi(klon)
77      ! cufi-----input-R- resolution des mailles en x (m)      ! cufi-----input-R- resolution des mailles en x (m)
78      ! cvfi-----input-R- resolution des mailles en y (m)      ! cvfi-----input-R- resolution des mailles en y (m)
79        REAL d_t(klon, klev), d_q(klon, klev)
80      ! d_t------output-R- le changement pour "t"      ! d_t------output-R- le changement pour "t"
81      ! d_q------output-R- le changement pour "q"      ! d_q------output-R- le changement pour "q"
82      ! d_u------output-R- le changement pour "u"  
83      ! d_v------output-R- le changement pour "v"      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
84      ! d_ts-----output-R- le changement pour "ts"      ! changement pour "u" et "v"
85    
86        REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)
87      ! 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)
88      !                    (orientation positive vers le bas)      !                    (orientation positive vers le bas)
89      ! 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)
90      ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal      REAL dflux_t(klon), dflux_q(klon)
     ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal  
91      ! dflux_t derive du flux sensible      ! dflux_t derive du flux sensible
92      ! dflux_q derive du flux latent      ! dflux_q derive du flux latent
93      !IM "slab" ocean      !IM "slab" ocean
94        REAL flux_o(klon), flux_g(klon)
95        !IM "slab" ocean
96      ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')      ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')
97      ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')      ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')
98        REAL y_flux_o(klon), y_flux_g(klon)
99        REAL tslab(klon), ytslab(klon)
100      ! tslab-in/output-R temperature du slab ocean (en Kelvin)      ! tslab-in/output-R temperature du slab ocean (en Kelvin)
101      ! uniqmnt pour slab      ! uniqmnt pour slab
102        REAL seaice(klon), y_seaice(klon)
103      ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')      ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')
104      !cc      REAL y_fqcalving(klon), y_ffonte(klon)
105        REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
106      ! ffonte----Flux thermique utilise pour fondre la neige      ! ffonte----Flux thermique utilise pour fondre la neige
107      ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la      ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
108      !           hauteur de neige, en kg/m2/s      !           hauteur de neige, en kg/m2/s
     ! on rajoute en output yu1 et yv1 qui sont les vents dans  
     ! la premiere couche  
     ! ces 4 variables sont maintenant traites dans phytrac  
     ! itr--------input-I- nombre de traceurs  
     ! tr---------input-R- q. de traceurs  
     ! flux_surf--input-R- flux de traceurs a la surface  
     ! d_tr-------output-R tendance de traceurs  
     !IM cf. AM : PBL  
     ! trmb1-------deep_cape  
     ! trmb2--------inhibition  
     ! trmb3-------Point Omega  
     ! Cape(klon)-------Cape du thermique  
     ! 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)  
     REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)  
109      REAL run_off_lic_0(klon), y_run_off_lic_0(klon)      REAL run_off_lic_0(klon), y_run_off_lic_0(klon)
110    
111      REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)      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
113        ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
114      REAL rugmer(klon), agesno(klon, nbsrf)      REAL rugmer(klon), agesno(klon, nbsrf)
115      REAL, INTENT (IN) :: rugoro(klon)      REAL, INTENT(IN):: rugoro(klon)
116      REAL cdragh(klon), cdragm(klon)      REAL, INTENT(out):: cdragh(klon), cdragm(klon)
     ! jour de l'annee en cours                  
     INTEGER jour  
     REAL rmu0(klon) ! cosinus de l'angle solaire zenithal      
117      ! taux CO2 atmosphere                          ! taux CO2 atmosphere                    
118      REAL co2_ppm      REAL co2_ppm
119      LOGICAL, INTENT (IN) :: debut      LOGICAL, INTENT(IN):: debut
120      LOGICAL, INTENT (IN) :: lafin      LOGICAL, INTENT(IN):: lafin
121      LOGICAL ok_veget      LOGICAL ok_veget
122      CHARACTER (len=*), INTENT (IN) :: ocean      CHARACTER(len=*), INTENT(IN):: ocean
123      INTEGER npas, nexca      INTEGER npas, nexca
124    
     REAL pctsrf(klon, nbsrf)  
125      REAL ts(klon, nbsrf)      REAL ts(klon, nbsrf)
126        ! ts-------input-R- temperature du sol (en Kelvin)
127      REAL d_ts(klon, nbsrf)      REAL d_ts(klon, nbsrf)
128        ! d_ts-----output-R- le changement pour "ts"
129      REAL snow(klon, nbsrf)      REAL snow(klon, nbsrf)
130      REAL qsurf(klon, nbsrf)      REAL qsurf(klon, nbsrf)
131      REAL evap(klon, nbsrf)      REAL evap(klon, nbsrf)
# Line 155  contains Line 134  contains
134    
135      REAL fluxlat(klon, nbsrf)      REAL fluxlat(klon, nbsrf)
136    
137      REAL rain_f(klon), snow_f(klon)      REAL, intent(in):: rain_fall(klon), snow_f(klon)
138      REAL fder(klon)      REAL fder(klon)
139    
140      REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)      REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)
141      REAL rugos(klon, nbsrf)      REAL rugos(klon, nbsrf)
142      ! la nouvelle repartition des surfaces sortie de l'interface      ! rugos----input-R- longeur de rugosite (en m)
     REAL pctsrf_new(klon, nbsrf)  
143    
144      REAL zcoefh(klon, klev)      REAL zcoefh(klon, klev)
145      REAL zu1(klon)      REAL zu1(klon)
146      REAL zv1(klon)      REAL zv1(klon)
147    
148      !$$$ PB ajout pour soil      !$$$ PB ajout pour soil
149      LOGICAL, INTENT (IN) :: soil_model      LOGICAL, INTENT(IN):: soil_model
150      !IM ajout seuils cdrm, cdrh      !IM ajout seuils cdrm, cdrh
151      REAL cdmmax, cdhmax      REAL cdmmax, cdhmax
152    
# Line 179  contains Line 157  contains
157      REAL ytsoil(klon, nsoilmx)      REAL ytsoil(klon, nsoilmx)
158      REAL qsol(klon)      REAL qsol(klon)
159    
     EXTERNAL clqh, clvent, coefkz, calbeta, cltrac  
   
160      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
161      REAL yalb(klon)      REAL yalb(klon)
162      REAL yalblw(klon)      REAL yalblw(klon)
163      REAL yu1(klon), yv1(klon)      REAL yu1(klon), yv1(klon)
164        ! on rajoute en output yu1 et yv1 qui sont les vents dans
165        ! la premiere couche
166      REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)      REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)
167      REAL yrain_f(klon), ysnow_f(klon)      REAL yrain_f(klon), ysnow_f(klon)
168      REAL ysollw(klon), ysolsw(klon), ysollwdown(klon)      REAL ysollw(klon), ysolsw(klon), ysollwdown(klon)
# Line 199  contains Line 177  contains
177      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)
178      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)
179      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
180      REAL ycoefh(klon, klev), ycoefm(klon, klev)      REAL coefh(klon, klev), coefm(klon, klev)
181      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
182      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
183      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)
# Line 208  contains Line 186  contains
186      PARAMETER (ok_nonloc=.FALSE.)      PARAMETER (ok_nonloc=.FALSE.)
187      REAL ycoefm0(klon, klev), ycoefh0(klon, klev)      REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
188    
     !IM 081204 hcl_Anne ? BEG  
189      REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)      REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
190      REAL ykmm(klon, klev+1), ykmn(klon, klev+1)      REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
191      REAL ykmq(klon, klev+1)      REAL ykmq(klon, klev+1)
192      REAL yq2(klon, klev+1), q2(klon, klev+1, nbsrf)      REAL yq2(klon, klev+1), q2(klon, klev+1, nbsrf)
193      REAL q2diag(klon, klev+1)      REAL q2diag(klon, klev+1)
     !IM 081204 hcl_Anne ? END  
194    
195      REAL u1lay(klon), v1lay(klon)      REAL u1lay(klon), v1lay(klon)
196      REAL delp(klon, klev)      REAL delp(klon, klev)
# Line 230  contains Line 206  contains
206    
207      ! maf pour sorties IOISPL en cas de debugagage      ! maf pour sorties IOISPL en cas de debugagage
208    
209      CHARACTER (80) cldebug      CHARACTER(80) cldebug
210      SAVE cldebug      SAVE cldebug
211      CHARACTER (8) cl_surf(nbsrf)      CHARACTER(8) cl_surf(nbsrf)
212      SAVE cl_surf      SAVE cl_surf
213      INTEGER nhoridbg, nidbg      INTEGER nhoridbg, nidbg
214      SAVE nhoridbg, nidbg      SAVE nhoridbg, nidbg
# Line 243  contains Line 219  contains
219      LOGICAL first_appel      LOGICAL first_appel
220      SAVE first_appel      SAVE first_appel
221      DATA first_appel/ .TRUE./      DATA first_appel/ .TRUE./
222      LOGICAL :: debugindex = .FALSE.      LOGICAL:: debugindex = .FALSE.
223      INTEGER idayref      INTEGER idayref
224      REAL t2m(klon, nbsrf), q2m(klon, nbsrf)      REAL t2m(klon, nbsrf), q2m(klon, nbsrf)
225      REAL u10m(klon, nbsrf), v10m(klon, nbsrf)      REAL u10m(klon, nbsrf), v10m(klon, nbsrf)
# Line 260  contains Line 236  contains
236      !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds      !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds
237      ! physiq ce qui permet de sortir les grdeurs par sous surface)      ! physiq ce qui permet de sortir les grdeurs par sous surface)
238      REAL pblh(klon, nbsrf)      REAL pblh(klon, nbsrf)
239        ! pblh------- HCL
240      REAL plcl(klon, nbsrf)      REAL plcl(klon, nbsrf)
241      REAL capcl(klon, nbsrf)      REAL capcl(klon, nbsrf)
242      REAL oliqcl(klon, nbsrf)      REAL oliqcl(klon, nbsrf)
243      REAL cteicl(klon, nbsrf)      REAL cteicl(klon, nbsrf)
244      REAL pblt(klon, nbsrf)      REAL pblt(klon, nbsrf)
245        ! pblT------- T au nveau HCL
246      REAL therm(klon, nbsrf)      REAL therm(klon, nbsrf)
247      REAL trmb1(klon, nbsrf)      REAL trmb1(klon, nbsrf)
248        ! trmb1-------deep_cape
249      REAL trmb2(klon, nbsrf)      REAL trmb2(klon, nbsrf)
250        ! trmb2--------inhibition
251      REAL trmb3(klon, nbsrf)      REAL trmb3(klon, nbsrf)
252        ! trmb3-------Point Omega
253      REAL ypblh(klon)      REAL ypblh(klon)
254      REAL ylcl(klon)      REAL ylcl(klon)
255      REAL ycapcl(klon)      REAL ycapcl(klon)
# Line 279  contains Line 260  contains
260      REAL ytrmb1(klon)      REAL ytrmb1(klon)
261      REAL ytrmb2(klon)      REAL ytrmb2(klon)
262      REAL ytrmb3(klon)      REAL ytrmb3(klon)
     REAL y_cd_h(klon), y_cd_m(klon)  
263      REAL uzon(klon), vmer(klon)      REAL uzon(klon), vmer(klon)
264      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
265      REAL psfce(klon), patm(klon)      REAL psfce(klon), patm(klon)
# Line 295  contains Line 275  contains
275      REAL t_coup      REAL t_coup
276      PARAMETER (t_coup=273.15)      PARAMETER (t_coup=273.15)
277    
278      CHARACTER (len=20) :: modname = 'clmain'      CHARACTER(len=20):: modname = 'clmain'
279    
280      !------------------------------------------------------------      !------------------------------------------------------------
281    
# Line 411  contains Line 391  contains
391      pctsrf_pot(:, is_oce) = 1. - zmasq      pctsrf_pot(:, is_oce) = 1. - zmasq
392      pctsrf_pot(:, is_sic) = 1. - zmasq      pctsrf_pot(:, is_sic) = 1. - zmasq
393    
394      DO nsrf = 1, nbsrf      loop_surface: DO nsrf = 1, nbsrf
395         ! chercher les indices:         ! Chercher les indices :
396         ni = 0         ni = 0
397         knon = 0         knon = 0
398         DO i = 1, klon         DO i = 1, klon
# Line 436  contains Line 416  contains
416            CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)            CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)
417         END IF         END IF
418    
419         IF (knon==0) CYCLE         if_knon: IF (knon /= 0) then
   
        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  
420            DO j = 1, knon            DO j = 1, knon
421               i = ni(j)               i = ni(j)
422               ypaprs(j, k) = paprs(i, k)               ypct(j) = pctsrf(i, nsrf)
423               ypplay(j, k) = pplay(i, k)               yts(j) = ts(i, nsrf)
424               ydelp(j, k) = delp(i, k)               ytslab(i) = tslab(i)
425               yu(j, k) = u(i, k)               ysnow(j) = snow(i, nsrf)
426               yv(j, k) = v(i, k)               yqsurf(j) = qsurf(i, nsrf)
427               yt(j, k) = t(i, k)               yalb(j) = albe(i, nsrf)
428               yq(j, k) = q(i, k)               yalblw(j) = alblw(i, nsrf)
429                 yrain_f(j) = rain_fall(i)
430                 ysnow_f(j) = snow_f(i)
431                 yagesno(j) = agesno(i, nsrf)
432                 yfder(j) = fder(i)
433                 ytaux(j) = flux_u(i, 1, nsrf)
434                 ytauy(j) = flux_v(i, 1, nsrf)
435                 ysolsw(j) = solsw(i, nsrf)
436                 ysollw(j) = sollw(i, nsrf)
437                 ysollwdown(j) = sollwdown(i)
438                 yrugos(j) = rugos(i, nsrf)
439                 yrugoro(j) = rugoro(i)
440                 yu1(j) = u1lay(i)
441                 yv1(j) = v1lay(i)
442                 yrads(j) = ysolsw(j) + ysollw(j)
443                 ypaprs(j, klev+1) = paprs(i, klev+1)
444                 y_run_off_lic_0(j) = run_off_lic_0(i)
445                 yu10mx(j) = u10m(i, nsrf)
446                 yu10my(j) = v10m(i, nsrf)
447                 ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))
448            END DO            END DO
        END DO  
449    
450         ! calculer Cdrag et les coefficients d'echange            ! IF bucket model for continent, copy soil water content
451         CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,&            IF (nsrf == is_ter .AND. .NOT. ok_veget) THEN
452              yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)               DO j = 1, knon
453         !IM 081204 BEG                  i = ni(j)
454         !CR test                  yqsol(j) = qsol(i)
        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))  
455               END DO               END DO
456            END DO            ELSE
457         END IF               yqsol = 0.
458              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)  
459    
460            IF (1==1) THEN            DO k = 1, nsoilmx
461               DO k = 1, klev               DO j = 1, knon
462                  DO i = 1, knon                  i = ni(j)
463                     ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))                  ytsoil(j, k) = ftsoil(i, k, nsrf)
                    ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))  
                 END DO  
464               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  
465            END DO            END DO
466    
467            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  
468               DO j = 1, knon               DO j = 1, knon
469                  i = ni(j)                  i = ni(j)
470                  yq2(j, k) = q2(i, k, nsrf)                  ypaprs(j, k) = paprs(i, k)
471                    ypplay(j, k) = pplay(i, k)
472                    ydelp(j, k) = delp(i, k)
473                    yu(j, k) = u(i, k)
474                    yv(j, k) = v(i, k)
475                    yt(j, k) = t(i, k)
476                    yq(j, k) = q(i, k)
477               END DO               END DO
478            END DO            END DO
479    
480            !   Bug introduit volontairement pour converger avec les resultats            ! calculer Cdrag et les coefficients d'echange
481            !  du papier sur les thermiques.            CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &
482            IF (1==1) THEN                 yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))
483               y_cd_m(1:knon) = ycoefm(1:knon, 1)            IF (iflag_pbl == 1) THEN
484               y_cd_h(1:knon) = ycoefh(1:knon, 1)               CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
485            ELSE               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
486               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)  
487            END IF            END IF
           CALL ustarhb(knon, yu, yv, y_cd_m, yustar)  
488    
489            IF (prt_level>9) THEN            ! on seuille coefm et coefh
490               PRINT *, 'USTAR = ', yustar            IF (nsrf == is_oce) THEN
491                 coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
492                 coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
493            END IF            END IF
494    
495            !   iflag_pbl peut etre utilise comme longuer de melange            IF (ok_kzmin) THEN
496                 ! Calcul d'une diffusion minimale pour les conditions tres stables
497                 CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
498                      coefm(:, 1), ycoefm0, ycoefh0)
499                 coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
500                 coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
501               END IF
502    
503              IF (iflag_pbl >= 3) THEN
504                 ! MELLOR ET YAMADA adapté à Mars, Richard Fournier et
505                 ! Frédéric Hourdin
506                 yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
507                      + ypplay(:knon, 1))) &
508                      * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
509                 DO k = 2, klev
510                    yzlay(1:knon, k) = yzlay(1:knon, k-1) &
511                         + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
512                         / ypaprs(1:knon, k) &
513                         * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
514                 END DO
515                 DO k = 1, klev
516                    yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
517                         / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
518                 END DO
519                 yzlev(1:knon, 1) = 0.
520                 yzlev(:knon, klev+1) = 2. * yzlay(:knon, klev) &
521                      - yzlay(:knon, klev - 1)
522                 DO k = 2, klev
523                    yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
524                 END DO
525                 DO k = 1, klev + 1
526                    DO j = 1, knon
527                       i = ni(j)
528                       yq2(j, k) = q2(i, k, nsrf)
529                    END DO
530                 END DO
531    
532            IF (iflag_pbl>=11) THEN               CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
              CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &  
                   yu, yv, yteta, y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, &  
                   iflag_pbl)  
           ELSE  
              CALL yamada4(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, &  
                   yv, yteta, y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)  
           END IF  
533    
534            ycoefm(1:knon, 1) = y_cd_m(1:knon)               IF (prt_level > 9) THEN
535            ycoefh(1:knon, 1) = y_cd_h(1:knon)                  PRINT *, 'USTAR = ', yustar
536            ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)               END IF
537            ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)  
538         END IF               ! iflag_pbl peut être utilisé comme longueur de mélange
539    
540                 IF (iflag_pbl >= 11) THEN
541                    CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &
542                         yu, yv, yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, &
543                         yustar, iflag_pbl)
544                 ELSE
545                    CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
546                         coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
547                 END IF
548    
549         ! calculer la diffusion des vitesses "u" et "v"               coefm(:knon, 2:) = ykmm(:knon, 2:klev)
550         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &               coefh(:knon, 2:) = ykmn(:knon, 2:klev)
551              ydelp, y_d_u, y_flux_u)            END IF
552         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &  
553              ydelp, y_d_v, y_flux_v)            ! calculer la diffusion des vitesses "u" et "v"
554              CALL clvent(knon, dtime, yu1, yv1, coefm, yt, yu, ypaprs, ypplay, &
555         ! pour le couplage                 ydelp, y_d_u, y_flux_u)
556         ytaux = y_flux_u(:, 1)            CALL clvent(knon, dtime, yu1, yv1, coefm, yt, yv, ypaprs, ypplay, &
557         ytauy = y_flux_v(:, 1)                 ydelp, y_d_v, y_flux_v)
558    
559         ! calculer la diffusion de "q" et de "h"            ! pour le couplage
560         CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,&            ytaux = y_flux_u(:, 1)
561              cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil,&            ytauy = y_flux_v(:, 1)
562              yqsol, ok_veget, ocean, npas, nexca, rmu0, co2_ppm, yrugos,&  
563              yrugoro, yu1, yv1, ycoefh, yt, yq, yts, ypaprs, ypplay,&            ! calculer la diffusion de "q" et de "h"
564              ydelp, yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, &            CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat, &
565              yfder, ytaux, ytauy, ywindsp, ysollw, ysollwdown, ysolsw,&                 cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil, &
566              yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts,&                 yqsol, ok_veget, ocean, npas, nexca, rmu0, co2_ppm, yrugos, &
567              yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,&                 yrugoro, yu1, yv1, coefh, yt, yq, yts, ypaprs, ypplay, &
568              y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g,&                 ydelp, yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, &
569              ytslab, y_seaice)                 yfder, ytaux, ytauy, ywindsp, ysollw, ysollwdown, ysolsw, &
570                   yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts, &
571         ! calculer la longueur de rugosite sur ocean                 yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q, &
572         yrugm = 0.                 y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g, &
573         IF (nsrf==is_oce) THEN                 ytslab, y_seaice)
574    
575              ! calculer la longueur de rugosite sur ocean
576              yrugm = 0.
577              IF (nsrf == is_oce) THEN
578                 DO j = 1, knon
579                    yrugm(j) = 0.018*coefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
580                         0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))
581                    yrugm(j) = max(1.5E-05, yrugm(j))
582                 END DO
583              END IF
584            DO j = 1, knon            DO j = 1, knon
585               yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &               y_dflux_t(j) = y_dflux_t(j)*ypct(j)
586                    0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2))               y_dflux_q(j) = y_dflux_q(j)*ypct(j)
587               yrugm(j) = max(1.5E-05, yrugm(j))               yu1(j) = yu1(j)*ypct(j)
588                 yv1(j) = yv1(j)*ypct(j)
589            END DO            END DO
        END IF  
        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  
590    
591         DO k = 1, klev            DO k = 1, klev
592            DO j = 1, knon               DO j = 1, knon
593               i = ni(j)                  i = ni(j)
594               ycoefh(j, k) = ycoefh(j, k)*ypct(j)                  coefh(j, k) = coefh(j, k)*ypct(j)
595               ycoefm(j, k) = ycoefm(j, k)*ypct(j)                  coefm(j, k) = coefm(j, k)*ypct(j)
596               y_d_t(j, k) = y_d_t(j, k)*ypct(j)                  y_d_t(j, k) = y_d_t(j, k)*ypct(j)
597               y_d_q(j, k) = y_d_q(j, k)*ypct(j)                  y_d_q(j, k) = y_d_q(j, k)*ypct(j)
598               !§§§ PB                  flux_t(i, k, nsrf) = y_flux_t(j, k)
599               flux_t(i, k, nsrf) = y_flux_t(j, k)                  flux_q(i, k, nsrf) = y_flux_q(j, k)
600               flux_q(i, k, nsrf) = y_flux_q(j, k)                  flux_u(i, k, nsrf) = y_flux_u(j, k)
601               flux_u(i, k, nsrf) = y_flux_u(j, k)                  flux_v(i, k, nsrf) = y_flux_v(j, k)
602               flux_v(i, k, nsrf) = y_flux_v(j, k)                  y_d_u(j, k) = y_d_u(j, k)*ypct(j)
603               !$$$ 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)
604               !$$$ 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)  
605            END DO            END DO
        END DO  
606    
607         evap(:, nsrf) = -flux_q(:, 1, nsrf)            evap(:, nsrf) = -flux_q(:, 1, nsrf)
608    
609         albe(:, nsrf) = 0.            albe(:, nsrf) = 0.
610         alblw(:, nsrf) = 0.            alblw(:, nsrf) = 0.
611         snow(:, nsrf) = 0.            snow(:, nsrf) = 0.
612         qsurf(:, nsrf) = 0.            qsurf(:, nsrf) = 0.
613         rugos(:, nsrf) = 0.            rugos(:, nsrf) = 0.
614         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  
615            DO j = 1, knon            DO j = 1, knon
616               i = ni(j)               i = ni(j)
617               qsol(i) = yqsol(j)               d_ts(i, nsrf) = y_d_ts(j)
618                 albe(i, nsrf) = yalb(j)
619                 alblw(i, nsrf) = yalblw(j)
620                 snow(i, nsrf) = ysnow(j)
621                 qsurf(i, nsrf) = yqsurf(j)
622                 rugos(i, nsrf) = yz0_new(j)
623                 fluxlat(i, nsrf) = yfluxlat(j)
624                 IF (nsrf == is_oce) THEN
625                    rugmer(i) = yrugm(j)
626                    rugos(i, nsrf) = yrugm(j)
627                 END IF
628                 agesno(i, nsrf) = yagesno(j)
629                 fqcalving(i, nsrf) = y_fqcalving(j)
630                 ffonte(i, nsrf) = y_ffonte(j)
631                 cdragh(i) = cdragh(i) + coefh(j, 1)
632                 cdragm(i) = cdragm(i) + coefm(j, 1)
633                 dflux_t(i) = dflux_t(i) + y_dflux_t(j)
634                 dflux_q(i) = dflux_q(i) + y_dflux_q(j)
635                 zu1(i) = zu1(i) + yu1(j)
636                 zv1(i) = zv1(i) + yv1(j)
637            END DO            END DO
638         END IF            IF (nsrf == is_ter) THEN
639         IF (nsrf==is_lic) THEN               DO j = 1, knon
640                    i = ni(j)
641                    qsol(i) = yqsol(j)
642                 END DO
643              END IF
644              IF (nsrf == is_lic) THEN
645                 DO j = 1, knon
646                    i = ni(j)
647                    run_off_lic_0(i) = y_run_off_lic_0(j)
648                 END DO
649              END IF
650              !$$$ PB ajout pour soil
651              ftsoil(:, :, nsrf) = 0.
652              DO k = 1, nsoilmx
653                 DO j = 1, knon
654                    i = ni(j)
655                    ftsoil(i, k, nsrf) = ytsoil(j, k)
656                 END DO
657              END DO
658    
659            DO j = 1, knon            DO j = 1, knon
660               i = ni(j)               i = ni(j)
661               run_off_lic_0(i) = y_run_off_lic_0(j)               DO k = 1, klev
662                    d_t(i, k) = d_t(i, k) + y_d_t(j, k)
663                    d_q(i, k) = d_q(i, k) + y_d_q(j, k)
664                    d_u(i, k) = d_u(i, k) + y_d_u(j, k)
665                    d_v(i, k) = d_v(i, k) + y_d_v(j, k)
666                    zcoefh(i, k) = zcoefh(i, k) + coefh(j, k)
667                 END DO
668            END DO            END DO
669         END IF  
670         !$$$ PB ajout pour soil            !cc diagnostic t, q a 2m et u, v a 10m
671         ftsoil(:, :, nsrf) = 0.  
        DO k = 1, nsoilmx  
672            DO j = 1, knon            DO j = 1, knon
673               i = ni(j)               i = ni(j)
674               ftsoil(i, k, nsrf) = ytsoil(j, k)               uzon(j) = yu(j, 1) + y_d_u(j, 1)
675            END DO               vmer(j) = yv(j, 1) + y_d_v(j, 1)
676         END DO               tair1(j) = yt(j, 1) + y_d_t(j, 1)
677                 qair1(j) = yq(j, 1) + y_d_q(j, 1)
678                 zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
679                      1)))*(ypaprs(j, 1)-ypplay(j, 1))
680                 tairsol(j) = yts(j) + y_d_ts(j)
681                 rugo1(j) = yrugos(j)
682                 IF (nsrf == is_oce) THEN
683                    rugo1(j) = rugos(i, nsrf)
684                 END IF
685                 psfce(j) = ypaprs(j, 1)
686                 patm(j) = ypplay(j, 1)
687    
688         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)  
689            END DO            END DO
        END DO  
690    
691         !cc diagnostic t, q a 2m et u, v a 10m            CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &
692                   zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &
693                   yt10m, yq10m, yu10m, yustar)
694    
695         DO j = 1, knon            DO j = 1, knon
696            i = ni(j)               i = ni(j)
697            uzon(j) = yu(j, 1) + y_d_u(j, 1)               t2m(i, nsrf) = yt2m(j)
698            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  
699    
700         CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &               ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
701              tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &               u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
702              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)  
703    
704         END DO            END DO
705    
706         DO i = 1, knon            CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &
707            y_cd_h(i) = ycoefh(i, 1)                 y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &
708            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  
709    
        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  
710            DO j = 1, knon            DO j = 1, knon
              ! on projette sur la grille globale  
711               i = ni(j)               i = ni(j)
712               IF (pctsrf_new(i, is_oce)>epsfra) THEN               pblh(i, nsrf) = ypblh(j)
713                  flux_o(i) = y_flux_o(j)               plcl(i, nsrf) = ylcl(j)
714               ELSE               capcl(i, nsrf) = ycapcl(j)
715                  flux_o(i) = 0.               oliqcl(i, nsrf) = yoliqcl(j)
716               END IF               cteicl(i, nsrf) = ycteicl(j)
717                 pblt(i, nsrf) = ypblt(j)
718                 therm(i, nsrf) = ytherm(j)
719                 trmb1(i, nsrf) = ytrmb1(j)
720                 trmb2(i, nsrf) = ytrmb2(j)
721                 trmb3(i, nsrf) = ytrmb3(j)
722            END DO            END DO
        END IF  
723    
        IF (nsrf==is_sic) THEN  
724            DO j = 1, knon            DO j = 1, knon
725               i = ni(j)               DO k = 1, klev + 1
726               ! On pondère lorsque l'on fait le bilan au sol :                  i = ni(j)
727               ! flux_g(i) = y_flux_g(j)*ypct(j)                  q2(i, k, nsrf) = yq2(j, k)
728               IF (pctsrf_new(i, is_sic)>epsfra) THEN               END DO
                 flux_g(i) = y_flux_g(j)  
              ELSE  
                 flux_g(i) = 0.  
              END IF  
729            END DO            END DO
730              !IM "slab" ocean
731              IF (nsrf == is_oce) THEN
732                 DO j = 1, knon
733                    ! on projette sur la grille globale
734                    i = ni(j)
735                    IF (pctsrf_new(i, is_oce)>epsfra) THEN
736                       flux_o(i) = y_flux_o(j)
737                    ELSE
738                       flux_o(i) = 0.
739                    END IF
740                 END DO
741              END IF
742    
743              IF (nsrf == is_sic) THEN
744                 DO j = 1, knon
745                    i = ni(j)
746                    ! On pondère lorsque l'on fait le bilan au sol :
747                    IF (pctsrf_new(i, is_sic)>epsfra) THEN
748                       flux_g(i) = y_flux_g(j)
749                    ELSE
750                       flux_g(i) = 0.
751                    END IF
752                 END DO
753    
        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                                                        
754            END IF            END IF
755            !OCEAN                                                                  IF (ocean == 'slab  ') THEN
756         END IF               IF (nsrf == is_oce) THEN
757      END DO                  tslab(1:klon) = ytslab(1:klon)
758                    seaice(1:klon) = y_seaice(1:klon)
759                 END IF
760              END IF
761           end IF if_knon
762        END DO loop_surface
763    
764      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
     ! A rajouter: conservation de l'albedo  
765    
766      rugos(:, is_oce) = rugmer      rugos(:, is_oce) = rugmer
767      pctsrf = pctsrf_new      pctsrf = pctsrf_new

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

  ViewVC Help
Powered by ViewVC 1.1.21