/[lmdze]/trunk/phylmd/pbl_surface.f
ViewVC logotype

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

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

  ViewVC Help
Powered by ViewVC 1.1.21