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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21