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

Diff of /trunk/phylmd/Interface_surf/pbl_surface.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/phylmd/Interface_surf/pbl_surface.f revision 301 by guez, Thu Aug 2 17:23:07 2018 UTC
# Line 1  Line 1 
1  module clmain_m  module pbl_surface_m
2    
3    IMPLICIT NONE    IMPLICIT NONE
4    
5  contains  contains
6    
7    SUBROUTINE clmain(dtime, itap, date0, pctsrf, pctsrf_new, t, q, u, v,&    SUBROUTINE pbl_surface(pctsrf, t, q, u, v, julien, mu0, ftsol, cdmmax, &
8         jour, rmu0, co2_ppm, ok_veget, ocean, npas, nexca, ts,&         cdhmax, ftsoil, qsol, paprs, pplay, fsnow, qsurf, evap, falbe, fluxlat, &
9         soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil,&         rain_fall, snow_f, fsolsw, fsollw, frugs, agesno, rugoro, d_t, d_q, &
10         qsol, paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat,&         d_u, d_v, d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, &
11         rain_f, snow_f, solsw, sollw, sollwdown, fder, rlon, rlat, cufi,&         dflux_t, dflux_q, coefh, t2m, q2m, u10m_srf, v10m_srf, pblh, capcl, &
12         cvfi, rugos, debut, lafin, agesno, rugoro, d_t, d_q, d_u, d_v,&         oliqcl, cteicl, pblt, therm, plcl, fqcalving, ffonte, run_off_lic_0)
13         d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2,&  
14         dflux_t, dflux_q, zcoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh,&      ! From phylmd/clmain.F, version 1.6, 2005/11/16 14:47:19
15         capcl, oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl,&      ! Author: Z. X. Li (LMD/CNRS), date: 1993 Aug. 18th
16         fqcalving, ffonte, run_off_lic_0, flux_o, flux_g, tslab, seaice)      ! Objet : interface de couche limite (diffusion verticale)
17    
18      ! From phylmd/clmain.F, version 1.6 2005/11/16 14:47:19      ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
19        ! de la couche limite pour les traceurs se fait avec "cltrac" et
20      ! Tout ce qui a trait aux traceurs est dans phytrac maintenant.      ! ne tient pas compte de la diff\'erentiation des sous-fractions
21      ! Pour l'instant le calcul de la couche limite pour les traceurs      ! de sol.
22      ! se fait avec cltrac et ne tient pas compte de la différentiation  
23      ! des sous-fractions de sol.      use cdrag_m, only: cdrag
24        use clqh_m, only: clqh
25      ! Pour pouvoir extraire les coefficients d'échanges et le vent      use clvent_m, only: clvent
26      ! dans la première couche, trois champs supplémentaires ont été créés :      use coef_diff_turb_m, only: coef_diff_turb
27      ! zcoefh, zu1 et zv1. Pour l'instant nous avons moyenné les valeurs      USE conf_gcm_m, ONLY: lmt_pas
28      ! de ces trois champs sur les 4 sous-surfaces du modèle. Dans l'avenir      USE conf_phys_m, ONLY: iflag_pbl
29      ! si les informations des sous-surfaces doivent être prises en compte      USE dimphy, ONLY: klev, klon
30      ! il faudra sortir ces mêmes champs en leur ajoutant une dimension,      USE dimsoil, ONLY: nsoilmx
     ! c'est a dire nbsrf (nombre de sous-surfaces).  
   
     ! Auteur Z.X. Li (LMD/CNRS) date: 1993/08/18  
     ! Objet : interface de "couche limite" (diffusion verticale)  
   
     ! Arguments:  
     ! dtime----input-R- interval du temps (secondes)  
     ! itap-----input-I- numero du pas de temps  
     ! date0----input-R- jour initial  
     ! t--------input-R- temperature (K)  
     ! q--------input-R- vapeur d'eau (kg/kg)  
     ! u--------input-R- vitesse u  
     ! v--------input-R- vitesse v  
     ! ts-------input-R- temperature du sol (en Kelvin)  
     ! paprs----input-R- pression a intercouche (Pa)  
     ! pplay----input-R- pression au milieu de couche (Pa)  
     ! radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2  
     ! rlat-----input-R- latitude en degree  
     ! 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)  
   
     ! d_t------output-R- le changement pour "t"  
     ! d_q------output-R- le changement pour "q"  
     ! d_u------output-R- le changement pour "u"  
     ! d_v------output-R- le changement pour "v"  
     ! d_ts-----output-R- le changement pour "ts"  
     ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)  
     !                    (orientation positive vers le bas)  
     ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)  
     ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal  
     ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal  
     ! dflux_t derive du flux sensible  
     ! dflux_q derive du flux latent  
     !IM "slab" ocean  
     ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')  
     ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')  
   
     ! tslab-in/output-R temperature du slab ocean (en Kelvin)  
     ! uniqmnt pour slab  
   
     ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')  
     !cc  
     ! ffonte----Flux thermique utilise pour fondre la neige  
     ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la  
     !           hauteur de neige, en kg/m2/s  
     ! on rajoute en output yu1 et yv1 qui sont les vents dans  
     ! la premiere couche  
     ! ces 4 variables sont maintenant traites dans phytrac  
     ! itr--------input-I- nombre de traceurs  
     ! tr---------input-R- q. de traceurs  
     ! flux_surf--input-R- flux de traceurs a la surface  
     ! d_tr-------output-R tendance de traceurs  
     !IM cf. AM : PBL  
     ! trmb1-------deep_cape  
     ! trmb2--------inhibition  
     ! trmb3-------Point Omega  
     ! Cape(klon)-------Cape du thermique  
     ! EauLiq(klon)-------Eau liqu integr du thermique  
     ! ctei(klon)-------Critere d'instab d'entrainmt des nuages de CL  
     ! lcl------- Niveau de condensation  
     ! pblh------- HCL  
     ! pblT------- T au nveau HCL  
   
     USE histcom, ONLY : histbeg_totreg, histdef, histend, histsync  
     use histwrite_m, only: histwrite  
     use calendar, ONLY : ymds2ju  
     USE dimens_m, ONLY : iim, jjm  
     USE indicesol, ONLY : epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf  
     USE dimphy, ONLY : klev, klon, zmasq  
     USE dimsoil, ONLY : nsoilmx  
     USE temps, ONLY : annee_ref, itau_phy  
     USE dynetat0_m, ONLY : day_ini  
     USE iniprint, ONLY : prt_level  
     USE suphec_m, ONLY : rd, rg, rkappa  
     USE conf_phys_m, ONLY : iflag_pbl  
     USE gath_cpl, ONLY : gath2cpl  
31      use hbtm_m, only: hbtm      use hbtm_m, only: hbtm
32        USE histwrite_phy_m, ONLY: histwrite_phy
33      REAL, INTENT (IN) :: dtime      USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
34      REAL date0      USE interfoce_lim_m, ONLY: interfoce_lim
35      INTEGER, INTENT (IN) :: itap      use phyetat0_m, only: zmasq
36      REAL t(klon, klev), q(klon, klev)      use stdlevvar_m, only: stdlevvar
37      REAL u(klon, klev), v(klon, klev)      USE suphec_m, ONLY: rd, rg
38      REAL, INTENT (IN) :: paprs(klon, klev+1)      use time_phylmdz, only: itap
39      REAL, INTENT (IN) :: pplay(klon, klev)  
40      REAL, INTENT (IN) :: rlon(klon), rlat(klon)      REAL, INTENT(inout):: pctsrf(klon, nbsrf)
41      REAL cufi(klon), cvfi(klon)      ! tableau des pourcentages de surface de chaque maille
42      REAL d_t(klon, klev), d_q(klon, klev)  
43      REAL d_u(klon, klev), d_v(klon, klev)      REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
44      REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)      REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg / kg)
45      REAL dflux_t(klon), dflux_q(klon)      REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
46      !IM "slab" ocean      INTEGER, INTENT(IN):: julien ! jour de l'annee en cours
47      REAL flux_o(klon), flux_g(klon)      REAL, intent(in):: mu0(klon) ! cosinus de l'angle solaire zenithal    
48      REAL y_flux_o(klon), y_flux_g(klon)      REAL, INTENT(IN):: ftsol(:, :) ! (klon, nbsrf) temp\'erature du sol (en K)
49      REAL tslab(klon), ytslab(klon)      REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
50      REAL seaice(klon), y_seaice(klon)  
51      REAL y_fqcalving(klon), y_ffonte(klon)      REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
52      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)      ! soil temperature of surface fraction
53      REAL run_off_lic_0(klon), y_run_off_lic_0(klon)  
54        REAL, INTENT(inout):: qsol(:) ! (klon)
55      REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)      ! column-density of water in soil, in kg m-2
56      REAL rugmer(klon), agesno(klon, nbsrf)  
57      REAL, INTENT (IN) :: rugoro(klon)      REAL, INTENT(IN):: paprs(klon, klev + 1) ! pression a intercouche (Pa)
58      REAL cdragh(klon), cdragm(klon)      REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
59      ! jour de l'annee en cours                      REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse
     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)  
60      REAL qsurf(klon, nbsrf)      REAL qsurf(klon, nbsrf)
61      REAL evap(klon, nbsrf)      REAL evap(klon, nbsrf)
62      REAL albe(klon, nbsrf)      REAL, intent(inout):: falbe(klon, nbsrf)
63      REAL alblw(klon, nbsrf)      REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)
64    
65      REAL fluxlat(klon, nbsrf)      REAL, intent(in):: rain_fall(klon)
66        ! liquid water mass flux (kg / m2 / s), positive down
67    
68      REAL rain_f(klon), snow_f(klon)      REAL, intent(in):: snow_f(klon)
69      REAL fder(klon)      ! solid water mass flux (kg / m2 / s), positive down
70    
71      REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)      REAL, INTENT(IN):: fsolsw(klon, nbsrf), fsollw(klon, nbsrf)
72      REAL rugos(klon, nbsrf)      REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)
73      ! la nouvelle repartition des surfaces sortie de l'interface      real agesno(klon, nbsrf)
74      REAL pctsrf_new(klon, nbsrf)      REAL, INTENT(IN):: rugoro(klon)
75    
76      REAL zcoefh(klon, klev)      REAL, intent(out):: d_t(:, :), d_q(:, :) ! (klon, klev)
77      REAL zu1(klon)      ! changement pour t et q
     REAL zv1(klon)  
78    
79      !$$$ PB ajout pour soil      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
80      LOGICAL, INTENT (IN) :: soil_model      ! changement pour "u" et "v"
     !IM ajout seuils cdrm, cdrh  
     REAL cdmmax, cdhmax  
81    
82      REAL ksta, ksta_ter      REAL, intent(out):: d_ts(:, :) ! (klon, nbsrf) variation of ftsol
     LOGICAL ok_kzmin  
83    
84      REAL ftsoil(klon, nsoilmx, nbsrf)      REAL, intent(out):: flux_t(klon, nbsrf)
85      REAL ytsoil(klon, nsoilmx)      ! flux de chaleur sensible (Cp T) (W / m2) (orientation positive vers
86      REAL qsol(klon)      ! le bas) à la surface
87    
88        REAL, intent(out):: flux_q(klon, nbsrf)
89        ! flux de vapeur d'eau (kg / m2 / s) à la surface
90    
91        REAL, intent(out):: flux_u(klon, nbsrf), flux_v(klon, nbsrf)
92        ! tension du vent (flux turbulent de vent) à la surface, en Pa
93    
94        REAL, INTENT(out):: cdragh(klon), cdragm(klon)
95        real q2(klon, klev + 1, nbsrf)
96    
97        REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)
98        ! dflux_t derive du flux sensible
99        ! dflux_q derive du flux latent
100        ! IM "slab" ocean
101    
102        REAL, intent(out):: coefh(:, 2:) ! (klon, 2:klev)
103        ! Pour pouvoir extraire les coefficients d'\'echange, le champ
104        ! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
105        ! ce champ sur les quatre sous-surfaces du mod\`ele.
106    
107        REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)
108    
109        REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)
110        ! composantes du vent \`a 10m sans spirale d'Ekman
111    
112        ! Ionela Musat. Cf. Anne Mathieu : planetary boundary layer, hbtm.
113        ! Comme les autres diagnostics on cumule dans physiq ce qui permet
114        ! de sortir les grandeurs par sous-surface.
115        REAL pblh(klon, nbsrf) ! height of planetary boundary layer
116        REAL capcl(klon, nbsrf)
117        REAL oliqcl(klon, nbsrf)
118        REAL cteicl(klon, nbsrf)
119        REAL, INTENT(inout):: pblt(klon, nbsrf) ! T au nveau HCL
120        REAL therm(klon, nbsrf)
121        REAL plcl(klon, nbsrf)
122    
123      EXTERNAL clqh, clvent, coefkz, calbeta, cltrac      REAL, intent(out):: fqcalving(klon, nbsrf)
124        ! flux d'eau "perdue" par la surface et necessaire pour limiter la
125        ! hauteur de neige, en kg / m2 / s
126    
127      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      real ffonte(klon, nbsrf) ! flux thermique utilise pour fondre la neige
128        REAL, intent(inout):: run_off_lic_0(:) ! (klon)
129    
130        ! Local:
131    
132        ! la nouvelle repartition des surfaces sortie de l'interface
133        REAL, save:: pctsrf_new_oce(klon)
134        REAL, save:: pctsrf_new_sic(klon)
135    
136        REAL y_fqcalving(klon), y_ffonte(klon)
137        real y_run_off_lic_0(klon), y_run_off_lic(klon)
138        REAL run_off_lic(klon) ! ruissellement total
139        REAL rugmer(klon)
140        REAL ytsoil(klon, nsoilmx)
141        REAL yts(klon), ypct(klon), yz0_new(klon)
142        real yrugos(klon) ! longueur de rugosite (en m)
143      REAL yalb(klon)      REAL yalb(klon)
144      REAL yalblw(klon)      REAL snow(klon), yqsurf(klon), yagesno(klon)
145      REAL yu1(klon), yv1(klon)      real yqsol(klon) ! column-density of water in soil, in kg m-2
146      REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)      REAL yrain_f(klon) ! liquid water mass flux (kg / m2 / s), positive down
147      REAL yrain_f(klon), ysnow_f(klon)      REAL ysnow_f(klon) ! solid water mass flux (kg / m2 / s), positive down
     REAL ysollw(klon), ysolsw(klon), ysollwdown(klon)  
     REAL yfder(klon), ytaux(klon), ytauy(klon)  
148      REAL yrugm(klon), yrads(klon), yrugoro(klon)      REAL yrugm(klon), yrads(klon), yrugoro(klon)
   
149      REAL yfluxlat(klon)      REAL yfluxlat(klon)
   
150      REAL y_d_ts(klon)      REAL y_d_ts(klon)
151      REAL y_d_t(klon, klev), y_d_q(klon, klev)      REAL y_d_t(klon, klev), y_d_q(klon, klev)
152      REAL y_d_u(klon, klev), y_d_v(klon, klev)      REAL y_d_u(klon, klev), y_d_v(klon, klev)
153      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)      REAL y_flux_t(klon), y_flux_q(klon)
154      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)      REAL y_flux_u(klon), y_flux_v(klon)
155      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
156      REAL ycoefh(klon, klev), ycoefm(klon, klev)      REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
157        real ycdragh(klon), ycdragm(klon)
158      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
159      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
160      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)      REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
161        REAL yq2(klon, klev + 1)
     LOGICAL ok_nonloc  
     PARAMETER (ok_nonloc=.FALSE.)  
     REAL ycoefm0(klon, klev), ycoefh0(klon, klev)  
   
     !IM 081204 hcl_Anne ? BEG  
     REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)  
     REAL ykmm(klon, klev+1), ykmn(klon, klev+1)  
     REAL ykmq(klon, klev+1)  
     REAL yq2(klon, klev+1), q2(klon, klev+1, nbsrf)  
     REAL q2diag(klon, klev+1)  
     !IM 081204 hcl_Anne ? END  
   
     REAL u1lay(klon), v1lay(klon)  
162      REAL delp(klon, klev)      REAL delp(klon, klev)
163      INTEGER i, k, nsrf      INTEGER i, k, nsrf
   
164      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)  
165    
166      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.      REAL pctsrf_pot(klon, nbsrf)
167        ! "pourcentage potentiel" pour tenir compte des \'eventuelles
168      ! maf pour sorties IOISPL en cas de debugagage      ! apparitions ou disparitions de la glace de mer
169    
170      CHARACTER (80) cldebug      REAL yt2m(klon), yq2m(klon), wind10m(klon)
171      SAVE cldebug      REAL ustar(klon)
     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)  
   
     REAL yt2m(klon), yq2m(klon), yu10m(klon)  
     REAL yustar(klon)  
     ! -- LOOP  
     REAL yu10mx(klon)  
     REAL yu10my(klon)  
     REAL ywindsp(klon)  
     ! -- LOOP  
172    
173      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)  
174      REAL ypblh(klon)      REAL ypblh(klon)
175      REAL ylcl(klon)      REAL ylcl(klon)
176      REAL ycapcl(klon)      REAL ycapcl(klon)
# Line 275  contains Line 178  contains
178      REAL ycteicl(klon)      REAL ycteicl(klon)
179      REAL ypblt(klon)      REAL ypblt(klon)
180      REAL ytherm(klon)      REAL ytherm(klon)
181      REAL ytrmb1(klon)      REAL u1(klon), v1(klon)
     REAL ytrmb2(klon)  
     REAL ytrmb3(klon)  
     REAL y_cd_h(klon), y_cd_m(klon)  
     REAL uzon(klon), vmer(klon)  
182      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
183      REAL psfce(klon), patm(klon)      REAL psfce(klon), patm(klon)
184    
185      REAL qairsol(klon), zgeo1(klon)      REAL qairsol(klon), zgeo1(klon)
186      REAL rugo1(klon)      REAL rugo1(klon)
187        REAL zgeop(klon, klev)
     ! utiliser un jeu de fonctions simples                
     LOGICAL zxli  
     PARAMETER (zxli=.FALSE.)  
   
     REAL zt, zqs, zdelta, zcor  
     REAL t_coup  
     PARAMETER (t_coup=273.15)  
   
     CHARACTER (len=20) :: modname = 'clmain'  
188    
189      !------------------------------------------------------------      !------------------------------------------------------------
190    
     ! initialisation Anne  
191      ytherm = 0.      ytherm = 0.
192    
     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  
   
193      DO k = 1, klev ! epaisseur de couche      DO k = 1, klev ! epaisseur de couche
194         DO i = 1, klon         DO i = 1, klon
195            delp(i, k) = paprs(i, k) - paprs(i, k+1)            delp(i, k) = paprs(i, k) - paprs(i, k + 1)
196         END DO         END DO
197      END DO      END DO
     DO i = 1, klon ! vent de la premiere couche  
        zx_alf1 = 1.0  
        zx_alf2 = 1.0 - zx_alf1  
        u1lay(i) = u(i, 1)*zx_alf1 + u(i, 2)*zx_alf2  
        v1lay(i) = v(i, 1)*zx_alf1 + v(i, 2)*zx_alf2  
     END DO  
198    
199      ! initialisation:      ! Initialization:
200        rugmer = 0.
201      DO i = 1, klon      cdragh = 0.
202         rugmer(i) = 0.0      cdragm = 0.
203         cdragh(i) = 0.0      dflux_t = 0.
204         cdragm(i) = 0.0      dflux_q = 0.
205         dflux_t(i) = 0.0      ypct = 0.
206         dflux_q(i) = 0.0      yqsurf = 0.
207         zu1(i) = 0.0      yrain_f = 0.
208         zv1(i) = 0.0      ysnow_f = 0.
209      END DO      yrugos = 0.
210      ypct = 0.0      ypaprs = 0.
211      yts = 0.0      ypplay = 0.
212      ysnow = 0.0      ydelp = 0.
     yqsurf = 0.0  
     yalb = 0.0  
     yalblw = 0.0  
     yrain_f = 0.0  
     ysnow_f = 0.0  
     yfder = 0.0  
     ytaux = 0.0  
     ytauy = 0.0  
     ysolsw = 0.0  
     ysollw = 0.0  
     ysollwdown = 0.0  
     yrugos = 0.0  
     yu1 = 0.0  
     yv1 = 0.0  
     yrads = 0.0  
     ypaprs = 0.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  
     ytsoil = 999999.  
213      yrugoro = 0.      yrugoro = 0.
214      ! -- 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  
     yfluxlat = 0.  
215      flux_t = 0.      flux_t = 0.
216      flux_q = 0.      flux_q = 0.
217      flux_u = 0.      flux_u = 0.
218      flux_v = 0.      flux_v = 0.
219      DO k = 1, klev      fluxlat = 0.
220         DO i = 1, klon      d_t = 0.
221            d_t(i, k) = 0.0      d_q = 0.
222            d_q(i, k) = 0.0      d_u = 0.
223            d_u(i, k) = 0.0      d_v = 0.
224            d_v(i, k) = 0.0      coefh = 0.
225            zcoefh(i, k) = 0.0      fqcalving = 0.
226         END DO      run_off_lic = 0.
227      END DO  
228        ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
229      ! Boucler sur toutes les sous-fractions du sol:      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
230        ! (\`a affiner).
     ! Initialisation des "pourcentages potentiels". On considère ici qu'on  
     ! peut avoir potentiellement de la glace sur tout le domaine océanique  
     ! (à affiner)  
231    
232      pctsrf_pot = pctsrf      pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
233        pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
234      pctsrf_pot(:, is_oce) = 1. - zmasq      pctsrf_pot(:, is_oce) = 1. - zmasq
235      pctsrf_pot(:, is_sic) = 1. - zmasq      pctsrf_pot(:, is_sic) = 1. - zmasq
236    
237      DO nsrf = 1, nbsrf      ! Tester si c'est le moment de lire le fichier:
238         ! chercher les indices:      if (mod(itap - 1, lmt_pas) == 0) then
239           CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
240        endif
241    
242        ! Boucler sur toutes les sous-fractions du sol:
243    
244        loop_surface: DO nsrf = 1, nbsrf
245           ! Chercher les indices :
246         ni = 0         ni = 0
247         knon = 0         knon = 0
248         DO i = 1, klon         DO i = 1, klon
249            ! pour determiner le domaine a traiter on utilise les surfaces            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
250            ! "potentielles"            ! "potentielles"
251            IF (pctsrf_pot(i, nsrf) > epsfra) THEN            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
252               knon = knon + 1               knon = knon + 1
# Line 435  contains Line 254  contains
254            END IF            END IF
255         END DO         END DO
256    
257         ! 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  
258            DO j = 1, knon            DO j = 1, knon
259               i = ni(j)               i = ni(j)
260               ytsoil(j, k) = ftsoil(i, k, nsrf)               ypct(j) = pctsrf(i, nsrf)
261                 yts(j) = ftsol(i, nsrf)
262                 snow(j) = fsnow(i, nsrf)
263                 yqsurf(j) = qsurf(i, nsrf)
264                 yalb(j) = falbe(i, nsrf)
265                 yrain_f(j) = rain_fall(i)
266                 ysnow_f(j) = snow_f(i)
267                 yagesno(j) = agesno(i, nsrf)
268                 yrugos(j) = frugs(i, nsrf)
269                 yrugoro(j) = rugoro(i)
270                 yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)
271                 ypaprs(j, klev + 1) = paprs(i, klev + 1)
272                 y_run_off_lic_0(j) = run_off_lic_0(i)
273            END DO            END DO
        END DO  
        DO k = 1, klev  
           DO j = 1, knon  
              i = ni(j)  
              ypaprs(j, k) = paprs(i, k)  
              ypplay(j, k) = pplay(i, k)  
              ydelp(j, k) = delp(i, k)  
              yu(j, k) = u(i, k)  
              yv(j, k) = v(i, k)  
              yt(j, k) = t(i, k)  
              yq(j, k) = q(i, k)  
           END DO  
        END DO  
274    
275         ! calculer Cdrag et les coefficients d'echange            ! For continent, copy soil water content
276         CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,&            IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
             yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)  
        !IM 081204 BEG  
        !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))  
              END DO  
           END DO  
        END IF  
277    
278         !IM cf JLD : on seuille ycoefm et ycoefh            ytsoil(:knon, :) = ftsoil(ni(:knon), :, nsrf)
        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  
279    
        !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  
280            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  
281               DO j = 1, knon               DO j = 1, knon
282                  i = ni(j)                  i = ni(j)
283                  yq2(j, k) = q2(i, k, nsrf)                  ypaprs(j, k) = paprs(i, k)
284                    ypplay(j, k) = pplay(i, k)
285                    ydelp(j, k) = delp(i, k)
286                    yu(j, k) = u(i, k)
287                    yv(j, k) = v(i, k)
288                    yt(j, k) = t(i, k)
289                    yq(j, k) = q(i, k)
290               END DO               END DO
291            END DO            END DO
292    
293            !   Bug introduit volontairement pour converger avec les resultats            ! Calculer les géopotentiels de chaque couche:
           !  du papier sur les thermiques.  
           IF (1==1) THEN  
              y_cd_m(1:knon) = ycoefm(1:knon, 1)  
              y_cd_h(1:knon) = ycoefh(1:knon, 1)  
           ELSE  
              y_cd_h(1:knon) = ycoefm(1:knon, 1)  
              y_cd_m(1:knon) = ycoefh(1:knon, 1)  
           END IF  
           CALL ustarhb(knon, yu, yv, y_cd_m, yustar)  
294    
295            IF (prt_level>9) THEN            zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
296               PRINT *, 'USTAR = ', yustar                 + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
297    
298              DO k = 2, klev
299                 zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
300                      * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
301                      * (ypplay(:knon, k - 1) - ypplay(:knon, k))
302              ENDDO
303    
304              CALL cdrag(nsrf, sqrt(yu(:knon, 1)**2 + yv(:knon, 1)**2), &
305                   yt(:knon, 1), yq(:knon, 1), zgeop(:knon, 1), ypaprs(:knon, 1), &
306                   yts(:knon), yqsurf(:knon), yrugos(:knon), ycdragm(:knon), &
307                   ycdragh(:knon))
308    
309              IF (iflag_pbl == 1) THEN
310                 ycdragm(:knon) = max(ycdragm(:knon), 0.)
311                 ycdragh(:knon) = max(ycdragh(:knon), 0.)
312              end IF
313    
314              ! on met un seuil pour ycdragm et ycdragh
315              IF (nsrf == is_oce) THEN
316                 ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
317                 ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
318            END IF            END IF
319    
320            !   iflag_pbl peut etre utilise comme longuer de melange            IF (iflag_pbl >= 6) then
321                 DO k = 1, klev + 1
322                    DO j = 1, knon
323                       i = ni(j)
324                       yq2(j, k) = q2(i, k, nsrf)
325                    END DO
326                 END DO
327              end IF
328    
329              call coef_diff_turb(nsrf, ni(:knon), ypaprs(:knon, :), &
330                   ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
331                   yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
332                   ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
333    
334              CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
335                   ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
336                   ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
337                   y_flux_u(:knon))
338              CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
339                   ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
340                   ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
341                   y_flux_v(:knon))
342    
343              CALL clqh(julien, nsrf, ni(:knon), ytsoil(:knon, :), yqsol(:knon), &
344                   mu0(ni(:knon)), yrugos(:knon), yrugoro(:knon), yu(:knon, 1), &
345                   yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), yt(:knon, :), &
346                   yq(:knon, :), yts(:knon), ypaprs(:knon, :), ypplay(:knon, :), &
347                   ydelp(:knon, :), yrads(:knon), yalb(:knon), snow(:knon), &
348                   yqsurf(:knon), yrain_f(:knon), ysnow_f(:knon), yfluxlat(:knon), &
349                   pctsrf_new_sic(ni(:knon)), yagesno(:knon), y_d_t(:knon, :), &
350                   y_d_q(:knon, :), y_d_ts(:knon), yz0_new(:knon), &
351                   y_flux_t(:knon), y_flux_q(:knon), y_dflux_t(:knon), &
352                   y_dflux_q(:knon), y_fqcalving(:knon), y_ffonte(:knon), &
353                   y_run_off_lic_0(:knon), y_run_off_lic(:knon))
354    
355              ! calculer la longueur de rugosite sur ocean
356    
357              yrugm = 0.
358    
359            IF (iflag_pbl>=11) THEN            IF (nsrf == is_oce) THEN
360               CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &               DO j = 1, knon
361                    yu, yv, yteta, y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, &                  yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
362                    iflag_pbl)                       / rg + 0.11 * 14E-6 &
363            ELSE                       / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
364               CALL yamada4(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, &                  yrugm(j) = max(1.5E-05, yrugm(j))
365                    yv, yteta, y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)               END DO
366            END IF            END IF
367    
368            ycoefm(1:knon, 1) = y_cd_m(1:knon)            DO k = 1, klev
369            ycoefh(1:knon, 1) = y_cd_h(1:knon)               DO j = 1, knon
370            ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)                  i = ni(j)
371            ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)                  y_d_t(j, k) = y_d_t(j, k) * ypct(j)
372         END IF                  y_d_q(j, k) = y_d_q(j, k) * ypct(j)
373                    y_d_u(j, k) = y_d_u(j, k) * ypct(j)
374         ! calculer la diffusion des vitesses "u" et "v"                  y_d_v(j, k) = y_d_v(j, k) * ypct(j)
375         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &               END DO
             ydelp, y_d_u, y_flux_u)  
        CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &  
             ydelp, y_d_v, y_flux_v)  
   
        ! pour le couplage  
        ytaux = y_flux_u(:, 1)  
        ytauy = y_flux_v(:, 1)  
   
        ! FH modif sur le cdrag temperature  
        !$$$PB : déplace dans clcdrag  
        !$$$      do i=1, knon  
        !$$$         ycoefh(i, 1)=ycoefm(i, 1)*0.8  
        !$$$      enddo  
   
        ! calculer la diffusion de "q" et de "h"  
        CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,&  
             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  
           DO j = 1, knon  
              yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &  
                   0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2))  
              yrugm(j) = max(1.5E-05, yrugm(j))  
376            END DO            END DO
        END IF  
        DO j = 1, knon  
           y_dflux_t(j) = y_dflux_t(j)*ypct(j)  
           y_dflux_q(j) = y_dflux_q(j)*ypct(j)  
           yu1(j) = yu1(j)*ypct(j)  
           yv1(j) = yv1(j)*ypct(j)  
        END DO  
377    
378         DO k = 1, klev            flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
379              flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
380              flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
381              flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
382    
383              evap(:, nsrf) = -flux_q(:, nsrf)
384    
385              falbe(:, nsrf) = 0.
386              fsnow(:, nsrf) = 0.
387              qsurf(:, nsrf) = 0.
388              frugs(:, nsrf) = 0.
389            DO j = 1, knon            DO j = 1, knon
390               i = ni(j)               i = ni(j)
391               ycoefh(j, k) = ycoefh(j, k)*ypct(j)               d_ts(i, nsrf) = y_d_ts(j)
392               ycoefm(j, k) = ycoefm(j, k)*ypct(j)               falbe(i, nsrf) = yalb(j)
393               y_d_t(j, k) = y_d_t(j, k)*ypct(j)               fsnow(i, nsrf) = snow(j)
394               y_d_q(j, k) = y_d_q(j, k)*ypct(j)               qsurf(i, nsrf) = yqsurf(j)
395               !§§§ PB               frugs(i, nsrf) = yz0_new(j)
396               flux_t(i, k, nsrf) = y_flux_t(j, k)               fluxlat(i, nsrf) = yfluxlat(j)
397               flux_q(i, k, nsrf) = y_flux_q(j, k)               IF (nsrf == is_oce) THEN
398               flux_u(i, k, nsrf) = y_flux_u(j, k)                  rugmer(i) = yrugm(j)
399               flux_v(i, k, nsrf) = y_flux_v(j, k)                  frugs(i, nsrf) = yrugm(j)
400               !$$$ PB        y_flux_t(j, k) = y_flux_t(j, k) * ypct(j)               END IF
401               !$$$ PB        y_flux_q(j, k) = y_flux_q(j, k) * ypct(j)               agesno(i, nsrf) = yagesno(j)
402               y_d_u(j, k) = y_d_u(j, k)*ypct(j)               fqcalving(i, nsrf) = y_fqcalving(j)
403               y_d_v(j, k) = y_d_v(j, k)*ypct(j)               ffonte(i, nsrf) = y_ffonte(j)
404               !$$$ PB        y_flux_u(j, k) = y_flux_u(j, k) * ypct(j)               cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
405               !$$$ PB        y_flux_v(j, k) = y_flux_v(j, k) * ypct(j)               cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
406            END DO               dflux_t(i) = dflux_t(i) + y_dflux_t(j) * ypct(j)
407         END DO               dflux_q(i) = dflux_q(i) + y_dflux_q(j) * ypct(j)
408              END DO
409              IF (nsrf == is_ter) THEN
410                 qsol(ni(:knon)) = yqsol(:knon)
411              else IF (nsrf == is_lic) THEN
412                 DO j = 1, knon
413                    i = ni(j)
414                    run_off_lic_0(i) = y_run_off_lic_0(j)
415                    run_off_lic(i) = y_run_off_lic(j)
416                 END DO
417              END IF
418    
419         evap(:, nsrf) = -flux_q(:, 1, nsrf)            ftsoil(:, :, nsrf) = 0.
420              ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)
421    
        albe(:, nsrf) = 0.  
        alblw(:, nsrf) = 0.  
        snow(:, nsrf) = 0.  
        qsurf(:, nsrf) = 0.  
        rugos(:, 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  
           DO j = 1, knon  
              i = ni(j)  
              qsol(i) = yqsol(j)  
           END DO  
        END IF  
        IF (nsrf==is_lic) THEN  
422            DO j = 1, knon            DO j = 1, knon
423               i = ni(j)               i = ni(j)
424               run_off_lic_0(i) = y_run_off_lic_0(j)               DO k = 1, klev
425                    d_t(i, k) = d_t(i, k) + y_d_t(j, k)
426                    d_q(i, k) = d_q(i, k) + y_d_q(j, k)
427                    d_u(i, k) = d_u(i, k) + y_d_u(j, k)
428                    d_v(i, k) = d_v(i, k) + y_d_v(j, k)
429                 END DO
430            END DO            END DO
431         END IF  
432         !$$$ PB ajout pour soil            forall (k = 2:klev) coefh(ni(:knon), k) &
433         ftsoil(:, :, nsrf) = 0.                 = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypct(:knon)
434         DO k = 1, nsoilmx  
435              ! diagnostic t, q a 2m et u, v a 10m
436    
437            DO j = 1, knon            DO j = 1, knon
438               i = ni(j)               i = ni(j)
439               ftsoil(i, k, nsrf) = ytsoil(j, k)               u1(j) = yu(j, 1) + y_d_u(j, 1)
440            END DO               v1(j) = yv(j, 1) + y_d_v(j, 1)
441         END DO               tair1(j) = yt(j, 1) + y_d_t(j, 1)
442                 qair1(j) = yq(j, 1) + y_d_q(j, 1)
443                 zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
444                      1))) * (ypaprs(j, 1)-ypplay(j, 1))
445                 tairsol(j) = yts(j) + y_d_ts(j)
446                 rugo1(j) = yrugos(j)
447                 IF (nsrf == is_oce) THEN
448                    rugo1(j) = frugs(i, nsrf)
449                 END IF
450                 psfce(j) = ypaprs(j, 1)
451                 patm(j) = ypplay(j, 1)
452    
453         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)  
454            END DO            END DO
        END DO  
   
        !cc diagnostic t, q a 2m et u, v a 10m  
   
        DO j = 1, knon  
           i = ni(j)  
           uzon(j) = yu(j, 1) + y_d_u(j, 1)  
           vmer(j) = yv(j, 1) + y_d_v(j, 1)  
           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)  
455    
456            qairsol(j) = yqsurf(j)            CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), qair1, &
457         END DO                 zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, &
458                   yq10m, wind10m(:knon), ustar(:knon))
459    
460         CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &            DO j = 1, knon
461              tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &               i = ni(j)
462              yu10m, yustar)               t2m(i, nsrf) = yt2m(j)
463         !IM 081204 END               q2m(i, nsrf) = yq2m(j)
   
        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)  
464    
465         END DO               u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
466                      / sqrt(u1(j)**2 + v1(j)**2)
467                 v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
468                      / sqrt(u1(j)**2 + v1(j)**2)
469              END DO
470    
471         DO i = 1, knon            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
472            y_cd_h(i) = ycoefh(i, 1)                 y_flux_q(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
473            y_cd_m(i) = ycoefm(i, 1)                 yq(:knon, :), ypblh(:knon), ycapcl, yoliqcl, ycteicl, ypblt, &
474         END DO                 ytherm, ylcl)
        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  
475    
        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  
476            DO j = 1, knon            DO j = 1, knon
              ! on projette sur la grille globale  
477               i = ni(j)               i = ni(j)
478               IF (pctsrf_new(i, is_oce)>epsfra) THEN               pblh(i, nsrf) = ypblh(j)
479                  flux_o(i) = y_flux_o(j)               plcl(i, nsrf) = ylcl(j)
480               ELSE               capcl(i, nsrf) = ycapcl(j)
481                  flux_o(i) = 0.               oliqcl(i, nsrf) = yoliqcl(j)
482               END IF               cteicl(i, nsrf) = ycteicl(j)
483                 pblt(i, nsrf) = ypblt(j)
484                 therm(i, nsrf) = ytherm(j)
485            END DO            END DO
        END IF  
486    
        IF (nsrf==is_sic) THEN  
487            DO j = 1, knon            DO j = 1, knon
488               i = ni(j)               DO k = 1, klev + 1
489               ! On pondère lorsque l'on fait le bilan au sol :                  i = ni(j)
490               ! flux_g(i) = y_flux_g(j)*ypct(j)                  q2(i, k, nsrf) = yq2(j, k)
491               IF (pctsrf_new(i, is_sic)>epsfra) THEN               END DO
                 flux_g(i) = y_flux_g(j)  
              ELSE  
                 flux_g(i) = 0.  
              END IF  
492            END DO            END DO
493           else
494         END IF            fsnow(:, nsrf) = 0.
495         !nsrf.EQ.is_sic                                                     end IF if_knon
496         IF (ocean=='slab  ') THEN      END DO loop_surface
           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  
497    
498      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
499      ! A rajouter: conservation de l'albedo      frugs(:, is_oce) = rugmer
500        pctsrf(:, is_oce) = pctsrf_new_oce
501        pctsrf(:, is_sic) = pctsrf_new_sic
502    
503      rugos(:, is_oce) = rugmer      CALL histwrite_phy("run_off_lic", run_off_lic)
     pctsrf = pctsrf_new  
504    
505    END SUBROUTINE clmain    END SUBROUTINE pbl_surface
506    
507  end module clmain_m  end module pbl_surface_m

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

  ViewVC Help
Powered by ViewVC 1.1.21