Changeset 98


Ignore:
Timestamp:
Jul 5, 2000, 4:58:04 PM (24 years ago)
Author:
lmdzadmin
Message:

Interface avec les differentes surface, version de travail.LF

Location:
LMDZ.3.3/branches/rel-LF/libf/phylmd
Files:
16 edited

Legend:

Unmodified
Added
Removed
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/YOMCST.h

    r2 r98  
    1 C A1.0 Fundamental constants
     1! A1.0 Fundamental constants
    22      REAL RPI,RCLUM,RHPLA,RKBOL,RNAVO
    3 C A1.1 Astronomical constants
     3! A1.1 Astronomical constants
    44      REAL RDAY,REA,REPSM,RSIYEA,RSIDAY,ROMEGA
    5 c A1.1.bis Constantes concernant l'orbite de la Terre:
     5! A1.1.bis Constantes concernant l'orbite de la Terre:
    66      REAL R_ecc, R_peri, R_incl
    7 C A1.2 Geoide
     7! A1.2 Geoide
    88      REAL RA,RG,R1SA
    9 C A1.3 Radiation
     9! A1.3 Radiation
    1010      REAL RSIGMA,RI0
    11 C A1.4 Thermodynamic gas phase
     11! A1.4 Thermodynamic gas phase
    1212      REAL R,RMD,RMV,RD,RV,RCPD,RCPV,RCVD,RCVV
    1313      REAL RKAPPA,RETV
    14 C A1.5,6 Thermodynamic liquid,solid phases
     14! A1.5,6 Thermodynamic liquid,solid phases
    1515      REAL RCW,RCS
    16 C A1.7 Thermodynamic transition of phase
     16! A1.7 Thermodynamic transition of phase
    1717      REAL RLVTT,RLSTT,RLMLT,RTT,RATM
    18 C A1.8 Curve of saturation
     18! A1.8 Curve of saturation
    1919      REAL RESTT,RALPW,RBETW,RGAMW,RALPS,RBETS,RGAMS
    2020      REAL RALPD,RBETD,RGAMD
    21 C
     21!
    2222      COMMON/YOMCST/RPI   ,RCLUM ,RHPLA ,RKBOL ,RNAVO
    2323     S      ,RDAY  ,REA   ,REPSM ,RSIYEA,RSIDAY,ROMEGA
     
    3131     S      ,RESTT ,RALPW ,RBETW ,RGAMW ,RALPS ,RBETS ,RGAMS
    3232     S      ,RALPD ,RBETD ,RGAMD
    33 C    ------------------------------------------------------------------
     33!    ------------------------------------------------------------------
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/clim.h

    r2 r98  
    4242      PARAMETER ( CLIM_Void    = 0  )
    4343      PARAMETER ( CLIM_MaxMod  = 8 )
    44       PARAMETER ( CLIM_MaxPort = 16 )
     44      PARAMETER ( CLIM_MaxPort = 40 )
    4545      PARAMETER ( CLIM_MaxSegments = 160 )
    4646      PARAMETER ( CLIM_MaxLink = CLIM_MaxMod * CLIM_MaxPort )
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F

    r86 r98  
    11      SUBROUTINE clmain(dtime,pctsrf,t,q,u,v,
    2      .                  soil_model,ts,soilcap,soilflux,
    3      .                  paprs,pplay,radsol,snow,qsol,
    4      .                  xlat, rugos,
     2     .                  ok_veget,ts,
     3     .                  paprs,pplay,radsol,snow,qsol,evap,albe,
     4     .                  rain_f, snow_f, solsw, sollw,
     5     .                  rlon, rlat, rugos,
     6     .                  debut, lafin,
    57     .                  d_t,d_q,d_u,d_v,d_ts,
    68     .                  flux_t,flux_q,flux_u,flux_v,cdragh,cdragm,
     
    4143c beta-----input-R- coefficient de l'evaporation reelle (0 a 1)
    4244c dif_grnd-input-R- coeff. de diffusion (chaleur) vers le sol profond
    43 c xlat-----input-R- latitude en degree
     45c rlat-----input-R- latitude en degree
    4446c rugos----input-R- longeur de rugosite (en m)
    4547c
     
    7577      REAL u(klon,klev), v(klon,klev)
    7678      REAL paprs(klon,klev+1), pplay(klon,klev), radsol(klon)
    77       REAL xlat(klon)
     79      REAL rlon(klon), rlat(klon)
    7880      REAL d_t(klon, klev), d_q(klon, klev)
    7981      REAL d_u(klon, klev), d_v(klon, klev)
    80       REAL flux_t(klon,klev), flux_q(klon,klev)
     82      REAL flux_t(klon,klev, nbsrf), flux_q(klon,klev, nbsrf)
    8183      REAL dflux_t(klon), dflux_q(klon)
    82       REAL flux_u(klon,klev), flux_v(klon,klev)
     84      REAL flux_u(klon,klev, nbsrf), flux_v(klon,klev, nbsrf)
    8385      REAL rugmer(klon)
    8486      REAL cdragh(klon), cdragm(klon)
     87      LOGICAL debut, lafin, ok_veget
    8588cAA      INTEGER itr
    8689cAA      REAL tr(klon,klev,nbtr)
     
    9396      REAL snow(klon,nbsrf)
    9497      REAL qsol(klon,nbsrf)
     98      REAL evap(klon,nbsrf)
     99      REAL albe(klon,nbsrf)
     100      real rain_f(klon), snow_f(klon)
     101      REAL sollw(klon), solsw(klon)
    95102      REAL rugos(klon,nbsrf)
    96103cAA
     
    103110c======================================================================
    104111      REAL yts(klon), yrugos(klon), ypct(klon)
    105       REAL ycal(klon), ybeta(klon), ydif(klon)
     112      REAL ycal(klon), ybeta(klon), ydif(klon), yalb(klon),yevap(klon)
    106113      REAL yu1(klon), yv1(klon)
     114      real ysnow(klon), yqsol(klon)
     115      real yrain_f(klon), ysnow_f(klon)
     116      real ysollw(klon), ysolsw(klon), ysolswnet(klon)
    107117      REAL yrugm(klon), yrads(klon)
    108118      REAL y_d_ts(klon)
     
    144154c======================================================================
    145155
    146       write(*,*)'CLMAIN.NEW'
    147 
    148156      DO k = 1, klev   ! epaisseur de couche
    149157      DO i = 1, klon
     
    174182         d_ts(i,nsrf) = 0.0
    175183      ENDDO
    176       ENDDO
     184      END DO
     185C§§§ PB
     186      flux_t = 0.
     187      flux_q = 0.
     188      flux_u = 0.
     189      flux_v = 0.
    177190      DO k = 1, klev
    178191      DO i = 1, klon
    179192         d_t(i,k) = 0.0
    180193         d_q(i,k) = 0.0
    181          flux_t(i,k) = 0.0
    182          flux_q(i,k) = 0.0
     194c$$$         flux_t(i,k) = 0.0
     195c$$$         flux_q(i,k) = 0.0
    183196         d_u(i,k) = 0.0
    184197         d_v(i,k) = 0.0
    185          flux_u(i,k) = 0.0
    186          flux_v(i,k) = 0.0
     198c$$$         flux_u(i,k) = 0.0
     199c$$$         flux_v(i,k) = 0.0
    187200         zcoefh(i,k) = 0.0
    188201      ENDDO
     
    203216c
    204217c prescrire les proprietes du sol:
    205       CALL calbeta(dtime,nsrf,snow,qsol, beta,capsol,dif_grnd)
    206       IF (.NOT.soil_model) THEN
    207          DO i = 1, klon
    208             cal(i) = RCPD * capsol(i)
    209             totalflu(i) = radsol(i)
    210          ENDDO
    211       ELSE
    212          DO i = 1, klon
    213             totalflu(i) = soilflux(i,nsrf) + radsol(i)
    214             IF (nsrf.EQ.is_oce) THEN
    215                cal(i) = 0.0
    216             ELSE
    217                cal(i) = RCPD / soilcap(i,nsrf)
    218             ENDIF
    219          ENDDO
    220       ENDIF
    221 c
     218c      CALL calbeta(dtime,nsrf,snow,qsol, beta,capsol,dif_grnd)
     219c      IF (.NOT.soil_model) THEN
     220c         DO i = 1, klon
     221c            cal(i) = RCPD * capsol(i)
     222c            totalflu(i) = radsol(i)
     223c         ENDDO
     224c      ELSE
     225c         DO i = 1, klon
     226c            totalflu(i) = soilflux(i,nsrf) + radsol(i)
     227c            IF (nsrf.EQ.is_oce) THEN
     228c               cal(i) = 0.0
     229c            ELSE
     230c               cal(i) = RCPD / soilcap(i,nsrf)
     231c            ENDIF
     232c         ENDDO
     233c      ENDIF
     234c
     235      totalflu = radsol
     236
    222237c chercher les indices:
    223238      DO j = 1, klon
     
    237252        ypct(j) = pctsrf(i,nsrf)
    238253        yts(j) = ts(i,nsrf)
     254        ysnow(j) = snow(i,nsrf)
     255        yevap(j) = evap(i,nsrf)
     256        yqsol(j) = qsol(i,nsrf)
     257        yalb(j) = albe(i,nsrf)
     258        yrain_f(j) = rain_f(i)
     259        ysnow_f(j) = snow_f(i)
     260        ysolsw(j) = solsw(i)
     261        ysollw(j) = sollw(i)
    239262        yrugos(j) = rugos(i,nsrf)
    240263        yu1(j) = u1lay(i)
    241264        yv1(j) = v1lay(i)
    242265        yrads(j) = totalflu(i)
    243         ycal(j) = cal(i)
    244         ybeta(j) = beta(i)
    245         ydif(j) = dif_grnd(i)
     266c        ycal(j) = cal(i)
     267c        ybeta(j) = beta(i)
     268c        ydif(j) = dif_grnd(i)
    246269        ypaprs(j,klev+1) = paprs(i,klev+1)
    247270      ENDDO
     
    259282      ENDDO
    260283c
    261 cAA      IF (itr.GE.1) THEN
    262 cAA      DO it = 1, itr
    263 cAA        DO k = 1, klev
    264 cAA        DO j = 1, knon
    265 cAA        i = ni(j)
    266 cAA          ytr(j,k,it) = tr(i,k,it)
    267 cAA        ENDDO
    268 cAA        ENDDO
    269 cAA        DO j = 1, knon
    270 cAA        i = ni(j)
    271 cAA          yflxsrf(j,it) = flux_surf(i,it)
    272 cAA        ENDDO
    273 cAA      ENDDO
    274 cAA      ENDIF
    275284c
    276285c calculer Cdrag et les coefficients d'echange
     
    287296      ENDDO
    288297c
    289 c parametrisation non-locale:
    290       IF (ok_nonloc) THEN
    291          DO i = 1, knon
    292             y_cd_h(i) = ycoefh(i,1)
    293             y_cd_m(i) = ycoefm(i,1)
    294          ENDDO
    295          CALL nonlocal(knon, ypaprs, ypplay,
    296      .        yts,ybeta,yu,yv,yt,yq,
    297      .        y_cd_h, y_cd_m, ycoefm0, ycoefh0, ygamt, ygamq)
    298          DO k = 1, klev
    299          DO i = 1, knon
    300             ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
    301             ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
    302          ENDDO
    303          ENDDO
    304       ELSE
    305          IF (.NOT.contreg) THEN
    306             DO k = 2, klev
    307             DO i = 1, knon
    308                ygamq(i,k) = 0.0
    309                ygamt(i,k) = 0.0
    310             ENDDO
    311             ENDDO
    312          ELSE
    313          DO k = 3, klev
    314          DO i = 1, knon
    315             ygamq(i,k) = 0.0
    316             ygamt(i,k) = -1.0E-03
    317          ENDDO
    318          ENDDO
    319          DO i = 1, knon
    320             ygamq(i,2) = 0.0
    321             ygamt(i,2) = -2.5E-03
    322          ENDDO
    323          ENDIF
    324       ENDIF
    325298c
    326299c calculer la diffusion de "q" et de "h"
    327       CALL clqh(knon, dtime, nsrf,yu1, yv1,
     300      CALL clqh(knon, dtime, nsrf, ni, pctsrf, rlon, rlat,
     301     e          yu1, yv1,
    328302     e          ycoefh,yt,yq,yts,ypaprs,ypplay,ydelp,yrads,
    329      e          ycal,ybeta,ydif,ygamt,ygamq,
     303     e          yevap,yalb, ysnow, yqsol, yrain_f, ysnow_f,
     304     e          ysollw, ysolsw,
    330305     s          y_d_t, y_d_q, y_d_ts,
    331306     s          y_flux_t, y_flux_q, y_dflux_t, y_dflux_q)
     
    363338c
    364339      DO k = 1, klev
     340        DO j = 1, knon
     341          i = ni(j)
     342          ycoefh(j,k) = ycoefh(j,k) * ypct(j)
     343          ycoefm(j,k) = ycoefm(j,k) * ypct(j)
     344          y_d_t(j,k) = y_d_t(j,k) * ypct(j)
     345          y_d_q(j,k) = y_d_q(j,k) * ypct(j)
     346C§§§ PB
     347          flux_t(i,k,nsrf) = y_flux_t(j,k)
     348          flux_q(i,k,nsrf) = y_flux_q(j,k)
     349          flux_u(i,k,nsrf) = y_flux_u(j,k)
     350          flux_v(i,k,nsrf) = y_flux_v(j,k)
     351c$$$ PB        y_flux_t(j,k) = y_flux_t(j,k) * ypct(j)
     352c$$$ PB        y_flux_q(j,k) = y_flux_q(j,k) * ypct(j)
     353          y_d_u(j,k) = y_d_u(j,k) * ypct(j)
     354          y_d_v(j,k) = y_d_v(j,k) * ypct(j)
     355c$$$ PB        y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)
     356c$$$ PB        y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)
     357        ENDDO
     358      ENDDO
     359
     360      evap(:,nsrf) = - flux_q(:,1,nsrf)
     361c
    365362      DO j = 1, knon
    366          ycoefh(j,k) = ycoefh(j,k) * ypct(j)
    367          ycoefm(j,k) = ycoefm(j,k) * ypct(j)
    368          y_d_t(j,k) = y_d_t(j,k) * ypct(j)
    369          y_d_q(j,k) = y_d_q(j,k) * ypct(j)
    370          y_flux_t(j,k) = y_flux_t(j,k) * ypct(j)
    371          y_flux_q(j,k) = y_flux_q(j,k) * ypct(j)
    372          y_d_u(j,k) = y_d_u(j,k) * ypct(j)
    373          y_d_v(j,k) = y_d_v(j,k) * ypct(j)
    374          y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)
    375          y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)
    376       ENDDO
    377       ENDDO
    378 c
    379       DO j = 1, knon
    380       i = ni(j)
     363         i = ni(j)
    381364         d_ts(i,nsrf) = y_d_ts(j)
     365         albe(i,nsrf) = yalb(j)
     366         snow(i,nsrf) = ysnow(j)
     367         qsol(i,nsrf) = yqsol(j)
    382368         rugmer(i) = yrugm(j)
    383369         cdragh(i) = cdragh(i) + ycoefh(j,1)
     
    400386         d_t(i,k) = d_t(i,k) + y_d_t(j,k)
    401387         d_q(i,k) = d_q(i,k) + y_d_q(j,k)
    402          flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
    403          flux_q(i,k) = flux_q(i,k) + y_flux_q(j,k)
     388c$$$ PB        flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
     389c$$$         flux_q(i,k) = flux_q(i,k) + y_flux_q(j,k)
    404390         d_u(i,k) = d_u(i,k) + y_d_u(j,k)
    405391         d_v(i,k) = d_v(i,k) + y_d_v(j,k)
    406          flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
    407          flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)
     392c$$$  PB       flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
     393c$$$         flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)
    408394         zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
    409395      ENDDO
     
    43041699999 CONTINUE
    431417c
     418
    432419      RETURN
    433420      END
    434       SUBROUTINE clqh(knon,dtime,nisurf,u1lay,v1lay,coef,
     421      SUBROUTINE clqh(knon,dtime,nisurf,knindex,pctsrf, rlon, rlat,
     422     e                u1lay,v1lay,coef,
    435423     e                t,q,ts,paprs,pplay,
    436      e                delp,radsol,cal,beta,dif_grnd, gamt,gamq,
     424     e                delp,radsol,evap,albedo,snow,qsol,
     425     e                precip_rain, precip_snow,
     426     e                lwdown, swdown,
    437427     s                d_t, d_q, d_ts, flux_t, flux_q,dflux_s,dflux_l)
    438428
     
    446436#include "dimensions.h"
    447437#include "dimphy.h"
     438#include "YOMCST.h"
     439#include "YOETHF.h"
     440#include "FCTTRE.h"
     441#include "indicesol.h"
    448442c Arguments:
    449443      INTEGER knon
     
    462456      REAL q(klon,klev)       ! humidite specifique (kg/kg)
    463457      REAL ts(klon)           ! temperature du sol (K)
     458      REAL evap(klon)         ! evaporation au sol
    464459      REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
    465460      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
    466461      REAL delp(klon,klev)    ! epaisseur de couche en pression (Pa)
    467462      REAL radsol(klon)       ! ray. net au sol (Solaire+IR) W/m2
     463      REAL albedo(klon)       ! albedo de la surface
     464      REAL snow(klon)         ! hauteur de neige
     465      REAL qsol(klon)         ! humidite de la surface
     466      real precip_rain(klon), precip_snow(klon)
     467      integer knindex(klon)
     468      real pctsrf(klon,nbsrf)
     469      real rlon(klon), rlat(klon)
    468470c
    469471      REAL d_t(klon,klev)     ! incrementation de "t"
     
    521523c======================================================================
    522524      REAL zcor, zdelta, zcvm5
    523 #include "YOMCST.h"
    524 #include "YOETHF.h"
    525 #include "FCTTRE.h"
    526 #include "indicesol.h"
     525      logical contreg
     526      parameter (contreg=.true.)
    527527c======================================================================
    528528c Rajout pour l'interface
     
    530530      integer jour
    531531      integer nisurf
    532       integer knindex(klon)
    533532      logical debut, lafin, ok_veget
    534       real rlon(klon), rlat(klon)
    535       real zlev(klon), zlflu(klon)
     533      real zlev1(klon)
    536534      real temp_air(klon), spechum(klon)
    537535      real hum_air(klon), ccanopy(klon)
    538536      real tq_cdrag(klon), petAcoef(klon), peqAcoef(klon)
    539537      real petBcoef(klon), peqBcoef(klon)
    540       real precip_rain(klon), precip_snow(klon)
    541538      real lwdown(klon), swnet(klon), swdown(klon), ps(klon)
    542539      real p1lay(klon)
     
    545542
    546543! Parametres de sortie
    547       real evap(klon), fluxsens(klon), fluxlat(klon)
     544      real fluxsens(klon), fluxlat(klon)
    548545      real tsol_rad(klon), tsurf_new(klon), alb_new(klon)
    549546      real emis_new(klon), z0_new(klon)
    550       real dflux_l(klon), dflux_s(klon)
    551547      real pctsrf_new(klon,nbsrf)
    552548c
     549
     550      if (.not. contreg) then
     551        do k = 2, klev
     552          do i = 1, knon
     553            gamq(i,k) = 0.0
     554            gamt(i,k) = 0.0
     555          enddo
     556        enddo
     557      else
     558        do k = 3, klev
     559          do i = 1, knon
     560            gamq(i,k)= 0.0
     561            gamt(i,k)=  -1.0e-03
     562          enddo
     563        enddo
     564        do i = 1, knon
     565          gamq(i,2) = 0.0
     566          gamt(i,2) = -2.5e-03
     567        enddo
     568      endif
    553569
    554570      DO i = 1, knon
     
    650666      spechum=q(:,1)
    651667      p1lay = pplay(:,1)
     668      zlev1 = delp(:,1)
     669      swnet = swdown * (1 - albedo)
     670c En attendant mieux
     671      hum_air = 0.
     672      ccanopy = 0.
     673      tq_cdrag = 0.
    652674
    653675      CALL interfsurf(itime, dtime, jour,
    654      . klon, nisurf, knon, knindex, rlon, rlat,
     676     . klon, iim, jjm, nisurf, knon, knindex, pctsrf, rlon, rlat,
    655677     . debut, lafin, ok_veget,
    656      . zlev, zlflu, u1lay, v1lay, temp_air, spechum, hum_air, ccanopy,
     678     . zlev1, u1lay, v1lay, temp_air, spechum, hum_air, ccanopy,
    657679     . tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef,
    658680     . precip_rain, precip_snow, lwdown, swnet, swdown,
    659      . ts, p1lay, cal, beta, coef1lay, psref, radsol, dif_grnd,
     681     . albedo, snow, qsol,
     682     . ts, p1lay, coef1lay, psref, radsol,
    660683     . ocean,
    661684     . evap, fluxsens, fluxlat, dflux_l, dflux_s,             
    662      . tsol_rad, tsurf_new, alb_new, emis_new, z0_new, pctsrf_new)
     685     . tsol_rad, tsurf_new, alb_new, emis_new, z0_new, pctsrf_new,
     686     . zmasq)
    663687
    664688      flux_t(:,1) = fluxsens
    665       flux_q(:,1) = fluxlat
     689      flux_q(:,1) = - evap
    666690      d_ts = tsurf_new - ts
     691      albedo = alb_new
    667692
    668693c==== une fois on a zx_h_ts, on peut faire l'iteration ========
     
    12281253      RETURN
    12291254      END
    1230       SUBROUTINE calbeta(dtime,indice,snow,qsol,
     1255      SUBROUTINE calbeta(dtime,indice,knon,snow,qsol,
    12311256     .                    vbeta,vcal,vdif)
    12321257      IMPLICIT none
     
    12381263c Calculer quelques parametres pour appliquer la couche limite
    12391264c ------------------------------------------------------------
    1240 #include "dimensions.h"
    1241 #include "dimphy.h"
     1265!#include "dimensions.h"
     1266!#include "dimphy.h"
    12421267#include "YOMCST.h"
    12431268#include "indicesol.h"
     
    12561281c
    12571282      REAL dtime
    1258       REAL snow(klon,nbsrf), qsol(klon,nbsrf)
    1259       INTEGER indice
    1260 C
    1261       REAL vbeta(klon)
    1262       REAL vcal(klon)
    1263       REAL vdif(klon)
     1283      REAL snow(knon,nbsrf), qsol(knon,nbsrf)
     1284      INTEGER indice, knon
     1285C
     1286      REAL vbeta(knon)
     1287      REAL vcal(knon)
     1288      REAL vdif(knon)
    12641289C
    12651290      IF (indice.EQ.is_oce) THEN
    1266       DO i = 1, klon
     1291      DO i = 1, knon
    12671292          vcal(i)   = 0.0
    12681293          vbeta(i)  = 1.0
     
    12721297c
    12731298      IF (indice.EQ.is_sic) THEN
    1274       DO i = 1, klon
     1299      DO i = 1, knon
    12751300          vcal(i) = calice
    12761301          IF (snow(i,is_sic) .GT. 0.0) vcal(i) = calsno
     
    12821307c
    12831308      IF (indice.EQ.is_ter) THEN
    1284       DO i = 1, klon
     1309      DO i = 1, knon
    12851310          vcal(i) = calsol
    12861311          IF (snow(i,is_ter) .GT. 0.0) vcal(i) = calsno
     
    12911316c
    12921317      IF (indice.EQ.is_lic) THEN
    1293       DO i = 1, klon
     1318      DO i = 1, knon
    12941319          vcal(i) = calice
    12951320          IF (snow(i,is_lic) .GT. 0.0) vcal(i) = calsno
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/condsurf.F

    r79 r98  
    3131
    3232      LOGICAL newlmt
    33       PARAMETER (newlmt=.FALSE.)
     33      PARAMETER (newlmt=.TRUE.)
    3434
    3535      INTEGER     nannemax
     
    114114      IF (newlmt) THEN
    115115c
    116 c Fraction "ocean":
    117       ierr = NF_INQ_VARID (nid, "FOCE", nvarid)
    118       IF (ierr .NE. NF_NOERR) THEN
    119          PRINT*, "condsurf: Le champ <FOCE> est absent"
    120          CALL abort
    121       ENDIF
    122 #ifdef NC_DOUBLE
    123       ierr = NF_GET_VARA_DOUBLE(nid,nvarid,debut,epais,pctsrf(1,is_oce))
    124 #else
    125       ierr = NF_GET_VARA_REAL(nid,nvarid,debut,epais,pctsrf(1,is_oce))
    126 #endif
    127       IF (ierr .NE. NF_NOERR) THEN
    128          PRINT*, "condsurf: Lecture echouee pour <FOCE>"
    129          CALL abort
    130       ENDIF
     116c$$$c Fraction "ocean":
     117c$$$      ierr = NF_INQ_VARID (nid, "FOCE", nvarid)
     118c$$$      IF (ierr .NE. NF_NOERR) THEN
     119c$$$         PRINT*, "condsurf: Le champ <FOCE> est absent"
     120c$$$         CALL abort
     121c$$$      ENDIF
     122c$$$#ifdef NC_DOUBLE
     123c$$$      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,debut,epais,pctsrf(1,is_oce))
     124c$$$#else
     125c$$$      ierr = NF_GET_VARA_REAL(nid,nvarid,debut,epais,pctsrf(1,is_oce))
     126c$$$#endif
     127c$$$      IF (ierr .NE. NF_NOERR) THEN
     128c$$$         PRINT*, "condsurf: Lecture echouee pour <FOCE>"
     129c$$$         CALL abort
     130c$$$      ENDIF
    131131c
    132132c Fraction "glace de mer":
     133c
     134c
    133135      ierr = NF_INQ_VARID (nid, "FSIC", nvarid)
    134136      IF (ierr .NE. NF_NOERR) THEN
     
    145147         CALL abort
    146148      ENDIF
    147 c
    148 c Fraction "terre":
    149       ierr = NF_INQ_VARID (nid, "FTER", nvarid)
    150       IF (ierr .NE. NF_NOERR) THEN
    151          PRINT*, "condsurf: Le champ <FTER> est absent"
    152          CALL abort
    153       ENDIF
    154 #ifdef NC_DOUBLE
    155       ierr = NF_GET_VARA_DOUBLE(nid,nvarid,debut,epais,pctsrf(1,is_ter))
    156 #else
    157       ierr = NF_GET_VARA_REAL(nid,nvarid,debut,epais,pctsrf(1,is_ter))
    158 #endif
    159       IF (ierr .NE. NF_NOERR) THEN
    160          PRINT*, "condsurf: Lecture echouee pour <FTER>"
    161          CALL abort
    162       ENDIF
    163 c
    164 c Fraction "glacier terre":
    165       ierr = NF_INQ_VARID (nid, "FLIC", nvarid)
    166       IF (ierr .NE. NF_NOERR) THEN
    167          PRINT*, "condsurf: Le champ <FLIC> est absent"
    168          CALL abort
    169       ENDIF
    170 #ifdef NC_DOUBLE
    171       ierr = NF_GET_VARA_DOUBLE(nid,nvarid,debut,epais,pctsrf(1,is_lic))
    172 #else
    173       ierr = NF_GET_VARA_REAL(nid,nvarid,debut,epais,pctsrf(1,is_lic))
    174 #endif
    175       IF (ierr .NE. 0) THEN
    176          PRINT*, "condsurf: Lecture echouee pour <FLIC>"
    177          CALL abort
    178       ENDIF
     149C
     150C positionnement % ocean libre et verification qu'il y a compatibilite des soussurfaces
     151      pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon))
     152     $    - pctsrf(1 : klon, is_sic)
     153      DO i = 1, klon
     154        IF ( pctsrf(i, is_sic) .GT. (1. - zmasq(i)) ) THEN
     155            WRITE(*,*) 'condsurf : sea-ice et masque pb en ', i,
     156     $     pctsrf(i, is_sic), (1. - zmasq(i))     
     157            pctsrf(i, is_sic) = (1. - zmasq(i))
     158            pctsrf(i, is_oce) = 0.
     159        ENDIF         
     160        IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
     161     $       pctsrf(i, is_oce) + pctsrf(i, is_sic)  - 1.) .GT. EPSFRA)
     162     $       THEN
     163             WRITE(*,*) 'physiq : pb sous surface au point ', i,
     164     $           pctsrf(i, 1 : nbsrf)
     165         ENDIF
     166      END DO
     167c
     168c$$$c Fraction "terre":
     169c$$$      ierr = NF_INQ_VARID (nid, "FTER", nvarid)
     170c$$$      IF (ierr .NE. NF_NOERR) THEN
     171c$$$         PRINT*, "condsurf: Le champ <FTER> est absent"
     172c$$$         CALL abort
     173c$$$      ENDIF
     174c$$$#ifdef NC_DOUBLE
     175c$$$      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,debut,epais,pctsrf(1,is_ter))
     176c$$$#else
     177c$$$      ierr = NF_GET_VARA_REAL(nid,nvarid,debut,epais,pctsrf(1,is_ter))
     178c$$$#endif
     179c$$$      IF (ierr .NE. NF_NOERR) THEN
     180c$$$         PRINT*, "condsurf: Lecture echouee pour <FTER>"
     181c$$$         CALL abort
     182c$$$      ENDIF
     183c$$$c
     184c$$$c Fraction "glacier terre":
     185c$$$      ierr = NF_INQ_VARID (nid, "FLIC", nvarid)
     186c$$$      IF (ierr .NE. NF_NOERR) THEN
     187c$$$         PRINT*, "condsurf: Le champ <FLIC> est absent"
     188c$$$         CALL abort
     189c$$$      ENDIF
     190c$$$#ifdef NC_DOUBLE
     191c$$$      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,debut,epais,pctsrf(1,is_lic))
     192c$$$#else
     193c$$$      ierr = NF_GET_VARA_REAL(nid,nvarid,debut,epais,pctsrf(1,is_lic))
     194c$$$#endif
     195c$$$      IF (ierr .NE. 0) THEN
     196c$$$         PRINT*, "condsurf: Lecture echouee pour <FLIC>"
     197c$$$         CALL abort
     198c$$$      ENDIF
    179199c
    180200      ELSE ! test sur newlmt
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/dimphy.h

    r83 r98  
    77      PARAMETER (nbtr=nqmx-2+1/(nqmx-1))
    88c-----------------------------------------------------------------------
     9      REAL zmasq(KLON)
     10      COMMON/terreoce/zmasq
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/inc_cpl.h

    r79 r98  
    11C
    2 C -- inc_cpl.h 
     2C -- inc_cpl.h   1998-04
    33C    **********
    44C@
    5 C@  Contents : variables describing field restart file names
     5C@  Contents : variables describing pipe and field names
    66C@  --------
    77C@
    8 C@ -- cl_write/cl_read  : for fields to write/READ
    9 C@ -- cl_f_write/cl_f_read  : for fields to write/read
     8C@ -- cl_write  : for fields to write
     9C@
     10C@ -- cl_read  : for fields to read
    1011C@
    1112C     -------------------------------------------------------------------
    1213C
     14      INTEGER jpread, jpwrit
     15      PARAMETER (jpread=0, jpwrit=1)
    1316      CHARACTER*8 cl_writ(jpmaxfld), cl_read(jpmaxfld)
    14       CHARACTER*8 cl_f_writ(jpmaxfld), cl_f_read(jpmaxfld)
     17      CHARACTER*6 cl_f_writ(jpmaxfld), cl_f_read(jpmaxfld)
    1518      COMMON / comcpl / cl_writ, cl_read, cl_f_writ, cl_f_read
    1619C     -------------------------------------------------------------------
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/inc_sipc.h

    r79 r98  
    1818      INTEGER  mpoolinitr
    1919      INTEGER  mpoolinitw
    20       INTEGER  mpoolwrit(jpflda2o)
    21       INTEGER  mpoolread(jpfldo2a)
     20      INTEGER  mpoolwrit(jpmaxfld)
     21      INTEGER  mpoolread(jpmaxfld)
    2222      COMMON / compool / mpoolinitr, mpoolinitw, mpoolwrit, mpoolread
    2323C     -------------------------------------------------------------------
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/indicesol.h

    r95 r98  
    1313      REAL epsfra
    1414      PARAMETER (epsfra=1.0E-05)
     15!
     16      CHARACTER *3 clnsurf(nbsrf)
     17      DATA clnsurf/'ter', 'lic', 'oce', 'sic'/
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90

    r97 r98  
    4545!
    4646  SUBROUTINE interfsurf_hq(itime, dtime, jour, &
    47       & klon, nisurf, knon, knindex, rlon, rlat, &
     47      & klon, iim, jjm, nisurf, knon, knindex, pctsrf, rlon, rlat, &
    4848      & debut, lafin, ok_veget, &
    49       & zlev, zlflu, u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, &
     49      & zlev, u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, &
    5050      & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
    5151      & precip_rain, precip_snow, lwdown, swnet, swdown, &
    52       & tsurf, p1lay, cal, beta, coef1lay, ps, radsol, dif_grnd, &
     52      & albedo, snow, qsol, &
     53      & tsurf, p1lay, coef1lay, ps, radsol, &
    5354      & ocean, &
    5455      & evap, fluxsens, fluxlat, dflux_l, dflux_s, &             
    55       & tsol_rad, tsurf_new, alb_new, emis_new, z0_new, pctsrf_new)
     56      & tsol_rad, tsurf_new, alb_new, emis_new, z0_new, pctsrf_new, zmasq)
    5657
    5758
     
    6768!   itime        numero du pas de temps
    6869!   klon         nombre total de points de grille
     70!   iim, jjm     nbres de pts de grille
    6971!   dtime        pas de temps de la physique (en s)
    7072!   jour         jour dans l'annee en cours
     
    7274!   knon         nombre de points de la surface a traiter
    7375!   knindex      index des points de la surface a traiter
     76!   pctsrf       tableau des pourcentages de surface de chaque maille
    7477!   rlon         longitudes
    7578!   rlat         latitudes
     
    7982!                     (si false calcul simplifie des fluxs sur les continents)
    8083!   zlev         hauteur de la premiere couche
    81 !   zlflu        epaisseur de la premier couche
     84!   u1_lay       vitesse u 1ere couche
     85!   v1_lay       vitesse v 1ere couche
     86!   temp_air     temperature de l'air 1ere couche
     87!   spechum      humidite specifique 1ere couche
     88!   hum_air      humidite de l'air
     89!   ccanopy      concentration CO2 canopee
     90!   tq_cdrag     cdrag
     91!   petAcoef     coeff. A de la resolution de la CL pour t
     92!   peqAcoef     coeff. A de la resolution de la CL pour q
     93!   petBcoef     coeff. B de la resolution de la CL pour t
     94!   peqBcoef     coeff. B de la resolution de la CL pour q
     95!   precip_rain  precipitation liquide
     96!   precip_snow  precipitation solide
     97!   lwdown       flux IR entrant a la surface
     98!   swnet        flux solaire net
     99!   swdown       flux solaire entrant a la surface
     100!   albedo       albedo de la surface
     101!   tsurf        temperature de surface
     102!   p1lay        pression 1er niveau (milieu de couche)
     103!   coef1lay     coefficient d'echange
     104!   ps           pression au sol
     105!   radsol       rayonnement net aus sol (LW + SW)
     106!   ocean        type d'ocean utilise (force, slab, couple)
     107!
     108! output:
     109!   evap         evaporation totale
     110!   fluxsens     flux de chaleur sensible
     111!   fluxlat      flux de chaleur latente
     112!   tsol_rad     
     113!   tsurf_new    temperature au sol
     114!   alb_new      albedo
     115!   emis_new     emissivite
     116!   z0_new       surface roughness
     117!   pctsrf_new   nouvelle repartition des surfaces
     118!   zmasq        masque terre/ocean
     119
     120  include 'indicesol.h'
     121
     122! Parametres d'entree
     123  integer, intent(IN) :: itime
     124  integer, intent(IN) :: iim, jjm
     125  integer, intent(IN) :: klon
     126  real, intent(IN) :: dtime
     127  integer, intent(IN) :: jour
     128  integer, intent(IN) :: nisurf
     129  integer, intent(IN) :: knon
     130  integer, dimension(knon), intent(in) :: knindex
     131  real, dimension(klon,nbsrf), intent(IN) :: pctsrf
     132  logical, intent(IN) :: debut, lafin, ok_veget
     133  real, dimension(klon), intent(IN) :: rlon, rlat
     134  real, dimension(knon), intent(IN) :: zlev
     135  real, dimension(knon), intent(IN) :: u1_lay, v1_lay
     136  real, dimension(knon), intent(IN) :: temp_air, spechum
     137  real, dimension(knon), intent(IN) :: hum_air, ccanopy
     138  real, dimension(knon), intent(IN) :: tq_cdrag, petAcoef, peqAcoef
     139  real, dimension(knon), intent(IN) :: petBcoef, peqBcoef
     140  real, dimension(knon), intent(IN) :: precip_rain, precip_snow
     141  real, dimension(knon), intent(IN) :: lwdown, swnet, swdown, ps, albedo
     142  real, dimension(knon), intent(IN) :: tsurf, p1lay, coef1lay
     143  real, dimension(knon), intent(IN) :: radsol
     144  real, dimension(klon), intent(IN) :: zmasq
     145  character (len = 6)  :: ocean
     146  real, dimension(knon), intent(INOUT) :: evap, snow, qsol
     147
     148! Parametres de sortie
     149  real, dimension(knon), intent(OUT):: fluxsens, fluxlat
     150  real, dimension(knon), intent(OUT):: tsol_rad, tsurf_new, alb_new
     151  real, dimension(knon), intent(OUT):: emis_new, z0_new
     152  real, dimension(knon), intent(OUT):: dflux_l, dflux_s
     153  real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new
     154
     155! Local
     156  character (len = 20) :: modname = 'interfsurf_hq'
     157  character (len = 80) :: abort_message
     158  logical, save        :: first_call = .true.
     159  integer              :: error
     160  logical              :: check = .true.
     161  real, dimension(knon):: cal, beta, dif_grnd, capsol
     162  real, parameter      :: calice=1.0/(5.1444e+06*0.15), tau_gl=1./86400.*5.
     163  real, parameter      :: calsno=1./(2.3867e+06*.15)
     164
     165#include "YOMCST.inc"
     166
     167  if (check) write(*,*) 'Entree ', modname
     168!
     169! On doit commencer par appeler les schemas de surfaces continentales
     170! car l'ocean a besoin du ruissellement qui est y calcule
     171!
     172  if (first_call) then
     173    if (nisurf /= is_ter .and. klon > 1) then
     174      write(*,*)' *** Warning ***'
     175      write(*,*)' nisurf = ',nisurf,' /= is_ter = ',is_ter
     176      write(*,*)'or on doit commencer par les surfaces continentales'
     177      abort_message='voir ci-dessus'
     178      call abort_gcm(modname,abort_message,1)
     179    endif
     180    if (ocean /= 'slab  ' .and. ocean /= 'force ' .and. ocean /= 'couple') then
     181      write(*,*)' *** Warning ***'
     182      write(*,*)'Option couplage pour l''ocean = ', ocean
     183      abort_message='option pour l''ocean non valable'
     184      call abort_gcm(modname,abort_message,1)
     185    endif
     186  endif
     187  first_call = .false.
     188!
     189! Calcul age de la neige
     190!
     191 
     192! Aiguillage vers les differents schemas de surface
     193
     194  if (nisurf == is_ter) then
     195!
     196! Surface "terre" appel a l'interface avec les sols continentaux
     197!
     198! allocation du run-off
     199    if (.not. allocated(run_off)) then
     200      allocate(run_off(knon), stat = error)
     201      if (error /= 0) then
     202        abort_message='Pb allocation run_off'
     203        call abort_gcm(modname,abort_message,1)
     204      endif
     205    else if (size(run_off) /= knon) then
     206      write(*,*)'Bizarre, le nombre de points continentaux'
     207      write(*,*)'a change entre deux appels. Je continue ...'
     208      deallocate(run_off, stat = error)
     209      allocate(run_off(knon), stat = error)
     210      if (error /= 0) then
     211        abort_message='Pb allocation run_off'
     212        call abort_gcm(modname,abort_message,1)
     213      endif
     214    endif
     215!
     216    if (.not. ok_veget) then
     217!
     218! calcul snow et qsol, hydrol adapté
     219!
     220      call calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
     221      cal = RCPD * capsol
     222      call calcul_fluxs( knon, nisurf, dtime, &
     223     &   tsurf, p1lay, cal, beta, coef1lay, ps, &
     224     &   precip_rain, precip_snow, snow, qsol,  &
     225     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
     226     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
     227     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
     228
     229!
     230! calcul albedo: lecture albedo fichier CL puis ajout albedo neige
     231!
     232    else
     233!
     234!  appel a sechiba
     235!
     236      call interfsol(itime, klon, dtime, nisurf, knon, &
     237     &  knindex, rlon, rlat, &
     238     &  debut, lafin, ok_veget, &
     239     &  zlev,  u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, &
     240     &  tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
     241     &  precip_rain, precip_snow, lwdown, swnet, swdown, &
     242     &  tsurf, p1lay, coef1lay, ps, radsol, &
     243     &  evap, fluxsens, fluxlat, &             
     244     &  tsol_rad, tsurf_new, alb_new, emis_new, z0_new, dflux_l, dflux_s)
     245    endif   
     246!
     247  else if (nisurf == is_oce) then
     248
     249    if (check) write(*,*)'ocean, nisurf = ',nisurf
     250
     251
     252!
     253! Surface "ocean" appel a l'interface avec l'ocean
     254!
     255!    if (ocean == 'couple') then
     256!      call interfoce(nisurf, ocean)
     257!    else if (ocean == 'slab  ') then
     258!      call interfoce(nisurf)
     259!    else                              ! lecture conditions limites
     260!      call interfoce(itime, dtime, jour, &
     261!     &  klon, nisurf, knon, knindex, &
     262!     &  debut, &
     263!     &  tsurf_new, alb_new, z0_new, pctsrf_new)
     264!    endif
     265!
     266    cal = 0.
     267    beta = 1.
     268    dif_grnd = 0.
     269
     270    call calcul_fluxs( knon, nisurf, dtime, &
     271     &   tsurf, p1lay, cal, beta, coef1lay, ps, &
     272     &   precip_rain, precip_snow, snow, qsol,  &
     273     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
     274     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
     275     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
     276
     277
     278!
     279! calcul albedo
     280!
     281
     282!
     283  else if (nisurf == is_sic) then
     284
     285    if (check) write(*,*)'sea ice, nisurf = ',nisurf
     286
     287!
     288! Surface "glace de mer" appel a l'interface avec l'ocean
     289!
     290!    call interfoce(nisurf, ocean)
     291!
     292
     293    cal = calice
     294    where (snow > 0.0) cal = calsno
     295    beta = 1.0
     296    dif_grnd = 1.0 / tau_gl
     297
     298    call calcul_fluxs( knon, nisurf, dtime, &
     299     &   tsurf, p1lay, cal, beta, coef1lay, ps, &
     300     &   precip_rain, precip_snow, snow, qsol,  &
     301     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
     302     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
     303     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
     304
     305
     306  else if (nisurf == is_lic) then
     307
     308    if (check) write(*,*)'glacier, nisurf = ',nisurf
     309
     310!
     311! Surface "glacier continentaux" appel a l'interface avec le sol
     312!
     313!    call interfsol(nisurf)
     314
     315    cal = calice
     316    where (snow > 0.0) cal = calsno
     317    beta = 1.0
     318    dif_grnd = 0.0
     319
     320    call calcul_fluxs( knon, nisurf, dtime, &
     321     &   tsurf, p1lay, cal, beta, coef1lay, ps, &
     322     &   precip_rain, precip_snow, snow, qsol,  &
     323     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
     324     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
     325     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
     326
     327  else
     328    write(*,*)'Index surface = ',nisurf
     329    abort_message = 'Index surface non valable'
     330    call abort_gcm(modname,abort_message,1)
     331  endif
     332
     333  END SUBROUTINE interfsurf_hq
     334
     335!
     336!#########################################################################
     337!
     338  SUBROUTINE interfsurf_vent(nisurf, knon   &         
     339  &                     )
     340!
     341! Cette routine sert d'aiguillage entre l'atmosphere et la surface en general
     342! (sols continentaux, oceans, glaces) pour les tensions de vents.
     343! En pratique l'interface se fait entre la couche limite du modele
     344! atmospherique (clmain.F) et les routines de surface (sechiba, oasis, ...)
     345!
     346!
     347! L.Fairhead 02/2000
     348!
     349! input:
     350!   nisurf       index de la surface a traiter (1 = sol continental)
     351!   knon         nombre de points de la surface a traiter
     352
     353! Parametres d'entree
     354  integer, intent(IN) :: nisurf
     355  integer, intent(IN) :: knon
     356
     357
     358  return
     359  END SUBROUTINE interfsurf_vent
     360!
     361!#########################################################################
     362!
     363  SUBROUTINE interfsol(itime, klon, dtime, nisurf, knon, &
     364     & knindex, rlon, rlat, &
     365     & debut, lafin, ok_veget, &
     366     & zlev,  u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, &
     367     & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
     368     & precip_rain, precip_snow, lwdown, swnet, swdown, &
     369     & tsurf, p1lay, coef1lay, ps, radsol, &
     370     & evap, fluxsens, fluxlat, &             
     371     & tsol_rad, tsurf_new, alb_new, emis_new, z0_new, dflux_l, dflux_s)
     372
     373! Cette routine sert d'interface entre le modele atmospherique et le
     374! modele de sol continental. Appel a sechiba
     375!
     376! L. Fairhead 02/2000
     377!
     378! input:
     379!   itime        numero du pas de temps
     380!   klon         nombre total de points de grille
     381!   dtime        pas de temps de la physique (en s)
     382!   nisurf       index de la surface a traiter (1 = sol continental)
     383!   knon         nombre de points de la surface a traiter
     384!   knindex      index des points de la surface a traiter
     385!   rlon         longitudes de la grille entiere
     386!   rlat         latitudes de la grille entiere
     387!   debut        logical: 1er appel a la physique (lire les restart)
     388!   lafin        logical: dernier appel a la physique (ecrire les restart)
     389!   ok_veget     logical: appel ou non au schema de surface continental
     390!                     (si false calcul simplifie des fluxs sur les continents)
     391!   zlev         hauteur de la premiere couche       
    82392!   u1_lay       vitesse u 1ere couche
    83393!   v1_lay       vitesse v 1ere couche
     
    98408!   tsurf        temperature de surface
    99409!   p1lay        pression 1er niveau (milieu de couche)
    100 !   cal          capacite calorifique du sol
    101 !   beta         evap reelle
    102410!   coef1lay     coefficient d'echange
    103411!   ps           pression au sol
    104412!   radsol       rayonnement net aus sol (LW + SW)
    105 !   dif_grnd     coeff. diffusion vers le sol profond
    106 !   ocean        type d'ocean utilise (force, slab, couple)
     413!   
     414!
     415! input/output
     416!   run_off      ruissellement total
    107417!
    108418! output:
     
    116426!   z0_new       surface roughness
    117427
    118   include 'indicesol.h'
    119 
    120 ! Parametres d'entree
    121   integer, intent(IN) :: itime
    122   integer, intent(IN) :: klon
    123   real, intent(IN) :: dtime
    124   integer, intent(IN) :: jour
    125   integer, intent(IN) :: nisurf
    126   integer, intent(IN) :: knon
    127   integer, dimension(knon), intent(in) :: knindex
    128   logical, intent(IN) :: debut, lafin, ok_veget
    129   real, dimension(klon), intent(IN) :: rlon, rlat
    130   real, dimension(knon), intent(IN) :: zlev, zlflu
    131   real, dimension(knon), intent(IN) :: u1_lay, v1_lay
    132   real, dimension(knon), intent(IN) :: temp_air, spechum
    133   real, dimension(knon), intent(IN) :: hum_air, ccanopy
    134   real, dimension(knon), intent(IN) :: tq_cdrag, petAcoef, peqAcoef
    135   real, dimension(knon), intent(IN) :: petBcoef, peqBcoef
    136   real, dimension(knon), intent(IN) :: precip_rain, precip_snow
    137   real, dimension(knon), intent(IN) :: lwdown, swnet, swdown, ps
    138   real, dimension(knon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
    139   real, dimension(knon), intent(IN) :: radsol, dif_grnd
    140   character (len = 6)  :: ocean
    141 
    142 ! Parametres de sortie
    143   real, dimension(knon), intent(OUT):: evap, fluxsens, fluxlat
    144   real, dimension(knon), intent(OUT):: tsol_rad, tsurf_new, alb_new
    145   real, dimension(knon), intent(OUT):: emis_new, z0_new
    146   real, dimension(knon), intent(OUT):: dflux_l, dflux_s
    147   real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new
    148 
    149 ! Local
    150   character (len = 20) :: modname = 'interfsurf_hq'
    151   character (len = 80) :: abort_message
    152   logical, save        :: first_call = .true.
    153   integer              :: error
    154   logical              :: check = .true.
    155 
    156   if (check) write(*,*) 'Entree ', modname
    157 !
    158 ! On doit commencer par appeler les schemas de surfaces continentales
    159 ! car l'ocean a besoin du ruissellement qui est y calcule
    160 !
    161   if (first_call) then
    162     if (nisurf /= is_ter .and. klon > 1) then
    163       write(*,*)' *** Warning ***'
    164       write(*,*)' nisurf = ',nisurf,' /= is_ter = ',is_ter
    165       write(*,*)'or on doit commencer par les surfaces continentales'
    166       abort_message='voir ci-dessus'
    167       call abort_gcm(modname,abort_message,1)
    168     endif
    169     if (ocean /= 'slab  ' .and. ocean /= 'force ' .and. ocean /= 'couple') then
    170       write(*,*)' *** Warning ***'
    171       write(*,*)'Option couplage pour l''ocean = ', ocean
    172       abort_message='option pour l''ocean non valable'
    173       call abort_gcm(modname,abort_message,1)
    174     endif
    175   endif
    176   first_call = .false.
    177 
    178 ! Aiguillage vers les differents schemas de surface
    179 
    180   if (nisurf == is_ter) then
    181 !
    182 ! Surface "terre" appel a l'interface avec les sols continentaux
    183 !
    184 ! allocation du run-off
    185     if (.not. allocated(run_off)) then
    186       allocate(run_off(knon), stat = error)
    187       if (error /= 0) then
    188         abort_message='Pb allocation run_off'
    189         call abort_gcm(modname,abort_message,1)
    190       endif
    191     else if (size(run_off) /= knon) then
    192       write(*,*)'Bizarre, le nombre de points continentaux'
    193       write(*,*)'a change entre deux appels. Je continue ...'
    194       deallocate(run_off, stat = error)
    195       allocate(run_off(knon), stat = error)
    196       if (error /= 0) then
    197         abort_message='Pb allocation run_off'
    198         call abort_gcm(modname,abort_message,1)
    199       endif
    200     endif
    201 !
    202     call interfsol(itime, klon, dtime, nisurf, knon, &
    203      &  knindex, rlon, rlat, &
    204      &  debut, lafin, ok_veget, &
    205      &  zlev, zlflu, u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, &
    206      &  tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
    207      &  precip_rain, precip_snow, lwdown, swnet, swdown, &
    208      &  tsurf, p1lay, cal, beta, coef1lay, ps, radsol, dif_grnd, &
    209      &  evap, fluxsens, fluxlat, &             
    210      &  tsol_rad, tsurf_new, alb_new, emis_new, z0_new, dflux_l, dflux_s)
    211      
    212 !
    213   else if (nisurf == is_oce) then
    214 
    215     if (check) write(*,*)'ocean, nisurf = ',nisurf
    216 !
    217 ! Surface "ocean" appel a l'interface avec l'ocean
    218 !
    219 !    if (ocean == 'couple') then
    220 !      call interfoce(nisurf, ocean)
    221 !    else if (ocean == 'slab  ') then
    222 !      call interfoce(nisurf)
    223 !    else                              ! lecture conditions limites
    224 !      call interfoce(itime, dtime, jour, &
    225 !     &  klon, nisurf, knon, knindex, &
    226 !     &  debut, &
    227 !     &  tsurf_new, alb_new, z0_new, pctsrf_new)
    228 !    endif
    229 !
    230     call calcul_fluxs( knon, dtime, &
    231      &   tsurf, p1lay, cal, beta, coef1lay, ps, &
    232      &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
    233      &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
    234      &   tsurf_new, fluxlat, fluxsens, dflux_s, dflux_l)
    235 !
    236   else if (nisurf == is_sic) then
    237 
    238     if (check) write(*,*)'sea ice, nisurf = ',nisurf
    239 
    240 !
    241 ! Surface "glace de mer" appel a l'interface avec l'ocean
    242 !
    243 !    call interfoce(nisurf, ocean)
    244 !
    245     call calcul_fluxs( knon, dtime, &
    246      &   tsurf, p1lay, cal, beta, coef1lay, ps, &
    247      &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
    248      &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
    249      &   tsurf_new, fluxlat, fluxsens, dflux_s, dflux_l)
    250 
    251   else if (nisurf == is_lic) then
    252 
    253     if (check) write(*,*)'glacier, nisurf = ',nisurf
    254 
    255 !
    256 ! Surface "glacier continentaux" appel a l'interface avec le sol
    257 !
    258 !    call interfsol(nisurf)
    259     call calcul_fluxs( knon, dtime, &
    260      &   tsurf, p1lay, cal, beta, coef1lay, ps, &
    261      &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
    262      &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
    263      &   tsurf_new, fluxlat, fluxsens, dflux_s, dflux_l)
    264 
    265   else
    266     write(*,*)'Index surface = ',nisurf
    267     abort_message = 'Index surface non valable'
    268     call abort_gcm(modname,abort_message,1)
    269   endif
    270 
    271   END SUBROUTINE interfsurf_hq
    272 
    273 !
    274 !#########################################################################
    275 !
    276   SUBROUTINE interfsurf_vent(nisurf, knon   &         
    277   &                     )
    278 !
    279 ! Cette routine sert d'aiguillage entre l'atmosphere et la surface en general
    280 ! (sols continentaux, oceans, glaces) pour les tensions de vents.
    281 ! En pratique l'interface se fait entre la couche limite du modele
    282 ! atmospherique (clmain.F) et les routines de surface (sechiba, oasis, ...)
    283 !
    284 !
    285 ! L.Fairhead 02/2000
    286 !
    287 ! input:
    288 !   nisurf       index de la surface a traiter (1 = sol continental)
    289 !   knon         nombre de points de la surface a traiter
    290 
    291 ! Parametres d'entree
    292   integer, intent(IN) :: nisurf
    293   integer, intent(IN) :: knon
    294 
    295 
    296   return
    297   END SUBROUTINE interfsurf_vent
    298 !
    299 !#########################################################################
    300 !
    301   SUBROUTINE interfsol(itime, klon, dtime, nisurf, knon, &
    302      & knindex, rlon, rlat, &
    303      & debut, lafin, ok_veget, &
    304      & zlev, zlflu, u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, &
    305      & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
    306      & precip_rain, precip_snow, lwdown, swnet, swdown, &
    307      & tsurf, p1lay, cal, beta, coef1lay, ps, radsol, dif_grnd, &
    308      & evap, fluxsens, fluxlat, &             
    309      & tsol_rad, tsurf_new, alb_new, emis_new, z0_new, dflux_l, dflux_s)
    310 
    311 ! Cette routine sert d'interface entre le modele atmospherique et le
    312 ! modele de sol continental. Appel a sechiba
    313 !
    314 ! L. Fairhead 02/2000
    315 !
    316 ! input:
    317 !   itime        numero du pas de temps
    318 !   klon         nombre total de points de grille
    319 !   dtime        pas de temps de la physique (en s)
    320 !   nisurf       index de la surface a traiter (1 = sol continental)
    321 !   knon         nombre de points de la surface a traiter
    322 !   knindex      index des points de la surface a traiter
    323 !   rlon         longitudes de la grille entiere
    324 !   rlat         latitudes de la grille entiere
    325 !   debut        logical: 1er appel a la physique (lire les restart)
    326 !   lafin        logical: dernier appel a la physique (ecrire les restart)
    327 !   ok_veget     logical: appel ou non au schema de surface continental
    328 !                     (si false calcul simplifie des fluxs sur les continents)
    329 !   zlev         hauteur de la premiere couche
    330 !   zlflu       
    331 !   u1_lay       vitesse u 1ere couche
    332 !   v1_lay       vitesse v 1ere couche
    333 !   temp_air     temperature de l'air 1ere couche
    334 !   spechum      humidite specifique 1ere couche
    335 !   hum_air      humidite de l'air
    336 !   ccanopy      concentration CO2 canopee
    337 !   tq_cdrag     cdrag
    338 !   petAcoef     coeff. A de la resolution de la CL pour t
    339 !   peqAcoef     coeff. A de la resolution de la CL pour q
    340 !   petBcoef     coeff. B de la resolution de la CL pour t
    341 !   peqBcoef     coeff. B de la resolution de la CL pour q
    342 !   precip_rain  precipitation liquide
    343 !   precip_snow  precipitation solide
    344 !   lwdown       flux IR entrant a la surface
    345 !   swnet        flux solaire net
    346 !   swdown       flux solaire entrant a la surface
    347 !   tsurf        temperature de surface
    348 !   p1lay        pression 1er niveau (milieu de couche)
    349 !   cal          capacite calorifique du sol
    350 !   beta         evap reelle
    351 !   coef1lay     coefficient d'echange
    352 !   ps           pression au sol
    353 !   radsol       rayonnement net aus sol (LW + SW)
    354 !   dif_grnd     coeff. diffusion vers le sol profond     
    355 !   
    356 !
    357 ! input/output
    358 !   run_off      ruissellement total
    359 !
    360 ! output:
    361 !   evap         evaporation totale
    362 !   fluxsens     flux de chaleur sensible
    363 !   fluxlat      flux de chaleur latente
    364 !   tsol_rad     
    365 !   tsurf_new    temperature au sol
    366 !   alb_new      albedo
    367 !   emis_new     emissivite
    368 !   z0_new       surface roughness
    369428
    370429! Parametres d'entree
     
    377436  logical, intent(IN) :: debut, lafin, ok_veget
    378437  real, dimension(klon), intent(IN) :: rlon, rlat
    379   real, dimension(knon), intent(IN) :: zlev, zlflu
     438  real, dimension(knon), intent(IN) :: zlev
    380439  real, dimension(knon), intent(IN) :: u1_lay, v1_lay
    381440  real, dimension(knon), intent(IN) :: temp_air, spechum
     
    385444  real, dimension(knon), intent(IN) :: precip_rain, precip_snow
    386445  real, dimension(knon), intent(IN) :: lwdown, swnet, swdown, ps
    387   real, dimension(knon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
    388   real, dimension(knon), intent(IN) :: radsol, dif_grnd
     446  real, dimension(knon), intent(IN) :: tsurf, p1lay, coef1lay
     447  real, dimension(knon), intent(IN) :: radsol
    389448! Parametres de sortie
    390449  real, dimension(knon), intent(OUT):: evap, fluxsens, fluxlat
     
    400459  character (len = 80) :: abort_message
    401460  logical              :: check = .true.
     461  real, dimension(knon) :: cal, beta, dif_grnd, capsol
    402462! type de couplage dans sechiba
    403463!  character (len=10)   :: coupling = 'implicit'
     
    414474  integer, save          :: rest_id_stom, hist_id_stom
    415475
     476  real, dimension(knon):: snow, qsol
     477
    416478  if (check) write(*,*)'Entree ', modname
    417479  if (check) write(*,*)'ok_veget = ',ok_veget
    418480
    419481! initialisation
    420   if (.not. ok_veget) then
    421     call calcul_fluxs( knon, dtime, &
    422      &   tsurf, p1lay, cal, beta, coef1lay, ps, &
    423      &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
    424      &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
    425      &   tsurf_new, fluxlat, fluxsens, dflux_s, dflux_l)
    426   else
    427482!    if (debut) then
    428483!      !
     
    476531!     & debut, lafin, coupling, control_in, &
    477532!     & lalo, neighbours, resolution,&
    478 !     & zlev, zlflu, u1_lay, v1_lay, spechum, temp_air,hum_air , ccanopy, &
     533!     & zlev, u1_lay, v1_lay, spechum, temp_air,hum_air , ccanopy, &
    479534!     & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
    480535!     & precip_rain, precip_snow, lwdown, swnet, swdown, ps, &
     
    486541! Sauvegarde dans fichiers histoire
    487542!
    488   endif  ! fin ok_veget
    489543
    490544  END SUBROUTINE interfsol
     
    492546!#########################################################################
    493547!
    494   SUBROUTINE interfoce_cpl(nisurf, ocean)
     548  SUBROUTINE interfoce_cpl(itime, dtime, &
     549      & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
     550      & ocean, nexca, debut, lafin, &
     551      & swdown, lwdown, precip_rain, precip_snow, evap, tsurf, &
     552      & fder, albsol, taux, tauy, &
     553      & tsurf_new, alb_new, alb_ice, pctsrf_new)
     554
     555
    495556
    496557! Cette routine sert d'interface entre le modele atmospherique et un
     
    500561!
    501562! input:
     563!   itime        numero du pas de temps
     564!   iim, jjm     nbres de pts de grille
     565!   dtime        pas de tempsde la physique
     566!   klon         nombre total de points de grille
    502567!   nisurf       index de la surface a traiter (1 = sol continental)
     568!   pctsrf       tableau des fractions de surface de chaque maille
     569!   knon         nombre de points de la surface a traiter
     570!   knindex      index des points de la surface a traiter
     571!   rlon         longitudes
     572!   rlat         latitudes
     573!   debut        logical: 1er appel a la physique
     574!   lafin        logical: dernier appel a la physique
    503575!   ocean        type d'ocean
     576!   nexca        frequence de couplage
     577!   swdown       flux solaire entrant a la surface
     578!   lwdown       flux IR entrant a la surface
     579!   precip_rain  precipitation liquide
     580!   precip_snow  precipitation solide
     581!   evap         evaporation
     582!   tsurf        temperature de surface
     583!   fder         derivee dF/dT
     584!   albsol       albedo du sol (coherent avec swdown)
     585!   taux         tension de vent en x
     586!   tauy         tension de vent en y
     587!   nexca        frequence de couplage
     588!
    504589!
    505590! output:
    506 !
     591!   tsurf_new    temperature au sol
     592!   alb_new      albedo
     593!   pctsrf_new   nouvelle repartition des surfaces
     594!   alb_ice      albedo de la glace
     595!
     596
     597#include 'indicesol.h'
    507598
    508599! Parametres d'entree
     600  integer, intent(IN) :: itime
     601  integer, intent(IN) :: iim, jjm
     602  real, intent(IN) :: dtime
     603  integer, intent(IN) :: klon
    509604  integer, intent(IN) :: nisurf
     605  integer, intent(IN) :: knon
     606  real, dimension(klon,nbsrf), intent(IN) :: pctsrf
     607  integer, dimension(knon), intent(in) :: knindex
     608  logical, intent(IN) :: debut, lafin
     609  real, dimension(klon), intent(IN) :: rlon, rlat
    510610  character (len = 6)  :: ocean
     611  real, dimension(knon), intent(IN) :: lwdown, swdown
     612  real, dimension(knon), intent(IN) :: precip_rain, precip_snow
     613  real, dimension(knon), intent(IN) :: tsurf, fder, albsol, taux, tauy
     614  integer              :: nexca
     615
     616  real, dimension(knon), intent(INOUT) :: evap
    511617
    512618! Parametres de sortie
     619  real, dimension(knon), intent(OUT):: tsurf_new, alb_new, alb_ice
     620  real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new
    513621
    514622! Variables locales
    515 
    516 
     623  integer                    :: j, error, sum_error, ig
     624  integer                    :: npas
     625  character (len = 20) :: modname = 'interfoce_cpl'
     626  character (len = 80) :: abort_message
     627  logical              :: check = .true.
     628! variables pour moyenner les variables de couplage
     629  real, allocatable, dimension(:),save :: cpl_sols, cpl_nsol, cpl_rain
     630  real, allocatable, dimension(:),save :: cpl_snow, cpl_evap, cpl_tsol
     631  real, allocatable, dimension(:),save :: cpl_fder, cpl_albe, cpl_taux
     632  real, allocatable, dimension(:),save :: cpl_tauy, cpl_ruis
     633! variables a passer au coupleur
     634  real, dimension(iim, jjm+1) :: wri_sols, wri_nsol, wri_rain
     635  real, dimension(iim, jjm+1) :: wri_snow, wri_evap, wri_tsol
     636  real, dimension(iim, jjm+1) :: wri_fder, wri_albe, wri_taux
     637  real, dimension(iim, jjm+1) :: wri_tauy, wri_ruis
     638! variables relues par le coupleur
     639  real, dimension(iim, jjm+1) :: read_sst, read_sic
     640  real, dimension(iim, jjm+1) :: read_alb_sst, read_alb_sic
     641! variable tampon
     642  real, dimension(klon)       :: tamp
     643  real, dimension(knon)       :: tamp_sic
     644
     645
     646!
    517647! Initialisation
     648!
     649  if (debut) then
     650    sum_error = 0
     651    allocate(cpl_sols(knon), stat = error)
     652    sum_error = sum_error + error
     653    allocate(cpl_nsol(knon), stat = error)
     654    sum_error = sum_error + error
     655    allocate(cpl_rain(knon), stat = error)
     656    sum_error = sum_error + error
     657    allocate(cpl_snow(knon), stat = error)
     658    sum_error = sum_error + error
     659    allocate(cpl_evap(knon), stat = error)
     660    sum_error = sum_error + error
     661    allocate(cpl_tsol(knon), stat = error)
     662    sum_error = sum_error + error
     663    allocate(cpl_fder(knon), stat = error)
     664    sum_error = sum_error + error
     665    allocate(cpl_albe(knon), stat = error)
     666    sum_error = sum_error + error
     667    allocate(cpl_taux(knon), stat = error)
     668    sum_error = sum_error + error
     669    allocate(cpl_tauy(knon), stat = error)
     670    sum_error = sum_error + error
     671    allocate(cpl_ruis(knon), stat = error)
     672    sum_error = sum_error + error
     673    if (sum_error /= 0) then
     674      abort_message='Pb allocation variables couplees'
     675      call abort_gcm(modname,abort_message,1)
     676    endif
     677    cpl_sols = 0.
     678    cpl_nsol = 0.
     679    cpl_rain = 0.
     680    cpl_snow = 0.
     681    cpl_evap = 0.
     682    cpl_tsol = 0.
     683    cpl_fder = 0.
     684    cpl_albe = 0.
     685    cpl_taux = 0.
     686    cpl_tauy = 0.
     687    cpl_ruis = 0.
     688!
     689! initialisation couplage
     690!
     691    call inicma(npas, nexca, dtime)
     692!
     693! 1ere lecture champs ocean
     694!
     695    call fromcpl(itime,(jjm+1)*iim,                                          &
     696     &        read_sst, read_sic, read_alb_sst, read_alb_sic)
     697    call cpl2gath(read_sst, tsurf_new, klon, knon,iim,jjm, knindex)
     698    call cpl2gath(read_sic, tamp_sic , klon, knon,iim,jjm, knindex)
     699    call cpl2gath(read_alb_sst, alb_new, klon, knon,iim,jjm, knindex)
     700    call cpl2gath(read_alb_sic, alb_ice, klon, knon,iim,jjm, knindex)
     701!
     702! transformer tamp_sic en pctsrf_new
     703!
     704    do ig = 1, klon
     705      IF (pctsrf(ig,is_oce) > epsfra .OR.            &
     706     &             pctsrf(ig,is_sic) > epsfra) THEN
     707            pctsrf_new(ig,is_oce) = pctsrf(ig,is_oce)    &
     708     &                        - (tamp_sic(ig)-pctsrf(ig,is_sic))
     709            pctsrf_new(ig,is_sic) = tamp_sic(ig)
     710      endif
     711    enddo
     712
     713  endif ! fin if (debut)
     714
    518715! fichier restart et fichiers histoires
    519716
    520717! calcul des fluxs a passer
    521718
     719  cpl_sols = cpl_sols + swdown      / FLOAT(nexca)
     720  cpl_nsol = cpl_nsol + lwdown      / FLOAT(nexca)
     721  cpl_rain = cpl_rain + precip_rain / FLOAT(nexca)
     722  cpl_snow = cpl_snow + precip_snow / FLOAT(nexca)
     723  cpl_evap = cpl_evap + evap        / FLOAT(nexca)
     724  cpl_tsol = cpl_tsol + tsurf       / FLOAT(nexca)
     725  cpl_fder = cpl_fder + fder        / FLOAT(nexca)
     726  cpl_albe = cpl_albe + albsol      / FLOAT(nexca)
     727  cpl_taux = cpl_taux + taux        / FLOAT(nexca)
     728  cpl_tauy = cpl_tauy + tauy        / FLOAT(nexca)
     729  cpl_ruis = cpl_ruis + run_off     / FLOAT(nexca)/dtime
     730
     731  if (mod(itime, nexca) == 0) then
     732!
     733! Mise sur la bonne grille des champs a passer au coupleur
     734    call gath2cpl(cpl_sols, wri_sols, klon, knon,iim,jjm, knindex)
     735    call gath2cpl(cpl_nsol, wri_nsol, klon, knon,iim,jjm, knindex)
     736    call gath2cpl(cpl_rain, wri_rain, klon, knon,iim,jjm, knindex)
     737    call gath2cpl(cpl_snow, wri_snow, klon, knon,iim,jjm, knindex)
     738    call gath2cpl(cpl_evap, wri_evap, klon, knon,iim,jjm, knindex)
     739    call gath2cpl(cpl_tsol, wri_tsol, klon, knon,iim,jjm, knindex)
     740    call gath2cpl(cpl_fder, wri_fder, klon, knon,iim,jjm, knindex)
     741    call gath2cpl(cpl_albe, wri_albe, klon, knon,iim,jjm, knindex)
     742    call gath2cpl(cpl_taux, wri_taux, klon, knon,iim,jjm, knindex)
     743    call gath2cpl(cpl_tauy, wri_tauy, klon, knon,iim,jjm, knindex)
     744    call gath2cpl(cpl_ruis, wri_ruis, klon, knon,iim,jjm, knindex)
     745!
     746! Passage des champs au coupleur
     747!
     748    call intocpl(itime, iim, jjm , wri_sols, wri_nsol, wri_rain, wri_snow, &
     749      & wri_evap, wri_tsol, wri_fder, wri_albe, wri_taux, wri_tauy,       &
     750      & wri_ruis )
     751    cpl_sols = 0.
     752    cpl_nsol = 0.
     753    cpl_rain = 0.
     754    cpl_snow = 0.
     755    cpl_evap = 0.
     756    cpl_tsol = 0.
     757    cpl_fder = 0.
     758    cpl_albe = 0.
     759    cpl_taux = 0.
     760    cpl_tauy = 0.
     761    cpl_ruis = 0.
     762
     763    call fromcpl(itime,(jjm+1)*iim,                                          &
     764     &        read_sst, read_sic, read_alb_sst, read_alb_sic)
     765    call cpl2gath(read_sst, tsurf_new, klon, knon,iim,jjm, knindex)
     766    call cpl2gath(read_sic, tamp_sic , klon, knon,iim,jjm, knindex)
     767    call cpl2gath(read_alb_sst, alb_new, klon, knon,iim,jjm, knindex)
     768    call cpl2gath(read_alb_sic, alb_ice, klon, knon,iim,jjm, knindex)
     769! transformer tamp_sic en pctsrf_new
     770
     771    do ig = 1, klon
     772      IF (pctsrf(ig,is_oce) > epsfra .OR.            &
     773     &             pctsrf(ig,is_sic) > epsfra) THEN
     774            pctsrf_new(ig,is_oce) = pctsrf(ig,is_oce)    &
     775     &                        - (tamp_sic(ig)-pctsrf(ig,is_sic))
     776            pctsrf_new(ig,is_sic) = tamp_sic(ig)
     777      endif
     778    enddo
     779  endif
     780
    522781  END SUBROUTINE interfoce_cpl
    523782!
    524783!#########################################################################
    525784!
     785
    526786  SUBROUTINE interfoce_slab(nisurf)
    527787
     
    8251085!
    8261086
    827   SUBROUTINE calcul_fluxs( knon, dtime, &
     1087  SUBROUTINE calcul_fluxs( knon, nisurf, dtime, &
    8281088     & tsurf, p1lay, cal, beta, coef1lay, ps, &
     1089     & precip_rain, precip_snow, snow, qsol, &
    8291090     & radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &
    8301091     & petAcoef, peqAcoef, petBcoef, peqBcoef, &
    831      & tsurf_new, fluxlat, fluxsens, dflux_s, dflux_l)
     1092     & tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
    8321093
    8331094! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
     
    8381099! input:
    8391100!   knon         nombre de points a traiter
     1101!   nisurf       surface a traiter
    8401102!   tsurf        temperature de surface
    8411103!   p1lay        pression 1er niveau (milieu de couche)
     
    8441106!   coef1lay     coefficient d'echange
    8451107!   ps           pression au sol
     1108!   precip_rain  precipitations liquides
     1109!   precip_snow  precipitations solides
     1110!   snow         champs hauteur de neige
     1111!   qsol         humidite du sol
     1112!   runoff       runoff en cas de trop plein
    8461113!   petAcoef     coeff. A de la resolution de la CL pour t
    8471114!   peqAcoef     coeff. A de la resolution de la CL pour q
     
    8621129#include "YOETHF.inc"
    8631130#include "FCTTRE.inc"
     1131#include 'indicesol.h'
    8641132
    8651133! Parametres d'entree
    866   integer, intent(IN) :: knon
     1134  integer, intent(IN) :: knon, nisurf
    8671135  real   , intent(IN) :: dtime
    8681136  real, dimension(knon), intent(IN) :: petAcoef, peqAcoef
     
    8701138  real, dimension(knon), intent(IN) :: ps, q1lay
    8711139  real, dimension(knon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
     1140  real, dimension(knon), intent(IN) :: precip_rain, precip_snow
    8721141  real, dimension(knon), intent(IN) :: radsol, dif_grnd
    8731142  real, dimension(knon), intent(IN) :: t1lay, u1lay, v1lay
     1143  real, dimension(knon), intent(INOUT) :: snow, qsol
    8741144
    8751145! Parametres sorties
    876   real, dimension(knon), intent(OUT):: tsurf_new, fluxlat, fluxsens
     1146  real, dimension(knon), intent(OUT):: tsurf_new, evap, fluxsens, fluxlat
    8771147  real, dimension(knon), intent(OUT):: dflux_s, dflux_l
    8781148
     
    8851155  real, dimension(knon) :: zx_h_ts, zx_q_0 , d_ts
    8861156  real                  :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
     1157  real                  :: bilan_f, fq_fonte
    8871158  real, parameter :: t_grnd = 271.35, t_coup = 273.15
    8881159  logical         :: check = .true.
    8891160  character (len = 20)  :: modname = 'calcul_fluxs'
     1161  logical         :: fonte_neige = .false.
     1162  real            :: max_eau_sol = 150.0
     1163  character (len = 80) :: abort_message
    8901164
    8911165  if (check) write(*,*)'Entree ', modname
     1166
     1167  if (size(run_off) /= knon .AND. nisurf == is_ter) then
     1168    write(*,*)'Bizarre, le nombre de points continentaux'
     1169    write(*,*)'a change entre deux appels. J''arrete ...'
     1170    abort_message='Pb run_off'
     1171    call abort_gcm(modname,abort_message,1)
     1172  endif
     1173!
     1174! Traitement neige et humidite du sol
     1175!
     1176    if (nisurf == is_oce) then
     1177      snow = 0.
     1178      qsol = max_eau_sol
     1179    else
     1180      snow = max(0.0, snow + (precip_snow - evap) * dtime)
     1181      qsol = qsol + (precip_rain - evap) * dtime
     1182    endif
     1183
     1184
    8921185!
    8931186! Initialisation
     
    9611254     &                       zx_nh(i) + zx_sl(i) * zx_nq(i)) & 
    9621255     &                     + dtime * dif_grnd(i))
     1256
     1257!
     1258! Y'a-t-il fonte de neige?
     1259!
     1260    fonte_neige = (nisurf /= is_oce) .AND. &
     1261     & (snow(i) > 0. .OR. nisurf == is_sic .OR. nisurf == is_lic) &
     1262     & .AND. (tsurf_new(i) >= RTT)
     1263    if (fonte_neige) tsurf_new(i) = RTT 
    9631264    zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
    9641265    d_ts(i) = tsurf_new(i) - tsurf(i)
     
    9661267!== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
    9671268!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
    968     fluxlat(i) = zx_mq(i) + zx_nq(i) * tsurf_new(i)
     1269    evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
     1270    fluxlat(i) = - evap(i) * zx_sl(i)
    9691271    fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
    9701272! Derives des flux dF/dTs (W m-2 K-1):
    9711273    dflux_s(i) = zx_nh(i)
    9721274    dflux_l(i) = (zx_sl(i) * zx_nq(i))
     1275!
     1276! en cas de fonte de neige
     1277!
     1278    if (fonte_neige) then
     1279      bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
     1280     &          dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
     1281     &          RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
     1282      bilan_f = max(0., bilan_f)
     1283      fq_fonte = bilan_f / zx_sl(i)
     1284      snow(i) = max(0., snow(i) - fq_fonte * dtime)
     1285      qsol(i) = qsol(i) + (fq_fonte * dtime)
     1286    endif
     1287    if (nisurf == is_ter)  &
     1288     &  run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
     1289    qsol(i) = min(qsol(i), max_eau_sol)
    9731290  ENDDO
    9741291
     
    9781295!
    9791296
     1297  SUBROUTINE sol_dem_write(itime, klon, rlon, rlat, &
     1298  &                     pctsrf_new,tsurf_new,alb_new)
     1299
     1300! Routine d'ecriture de l'etat de redemarrage pour le sol
     1301!
     1302! L.Fairhead
     1303!
     1304! input:
     1305!   itime        numero du pas de temps
     1306!   klon         nombre total de points de grille
     1307!   rlon         longitudes
     1308!   rlat         latitudes
     1309!   tsurf_new    temperature au sol
     1310!   alb_new      albedo
     1311!   pctsrf_new   repartition des surfaces
     1312
     1313  include 'indicesol.h'
     1314#include 'temps.inc'
     1315  include 'netcdf.inc'
     1316
     1317! Parametres d'entree
     1318  integer, intent(IN) :: itime
     1319  integer, intent(IN) :: klon
     1320  real, dimension(klon), intent(IN) :: rlon, rlat
     1321  real, dimension(klon,nbsrf), intent(IN)  :: tsurf_new, alb_new
     1322  real, dimension(klon,nbsrf), intent(IN) :: pctsrf_new
     1323
     1324! Variables locales
     1325  integer             :: ierr, nid
     1326  integer             :: idim1, idim2, idim3
     1327  integer,parameter   :: length = 100
     1328  character (len = 20) :: modname = 'sol_dem_write'
     1329  character (len = 80) :: abort_message
     1330  real, dimension(length) :: tab_cntrl = 0.
     1331  integer                 :: nvarid
     1332
     1333  ierr = NF_CREATE('restartsol', NF_CLOBBER, nid)
     1334  IF (ierr.NE.NF_NOERR) THEN
     1335    abort_message=' Pb d''ouverture du fichier restartsol'
     1336    CALL abort_gcm(modname,abort_message,ierr)
     1337  ENDIF
     1338
     1339  ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 23,    &
     1340     &                    "Fichier redemmarage sol")
     1341  ierr = NF_DEF_DIM (nid, "index", length, idim1)
     1342  ierr = NF_DEF_DIM (nid, "points_physiques", klon, idim2)
     1343  ierr = NF_DEF_DIM (nid, "nombre_surfaces", nbsrf, idim3)
     1344  ierr = NF_ENDDEF(nid)
     1345
     1346  tab_cntrl(13) = day_end
     1347  tab_cntrl(14) = anne_ini
     1348
     1349  ierr = NF_REDEF (nid)
     1350#ifdef NC_DOUBLE
     1351  ierr = NF_DEF_VAR (nid, "controle", NF_DOUBLE, 1, idim1,nvarid)
     1352#else
     1353  ierr = NF_DEF_VAR (nid, "controle", NF_FLOAT, 1, idim1,nvarid)
     1354#endif
     1355  ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 22,  &
     1356     &                        "Parametres de controle")
     1357  ierr = NF_ENDDEF(nid)
     1358#ifdef NC_DOUBLE
     1359  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
     1360#else
     1361  ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
     1362#endif
     1363
     1364  ierr = NF_REDEF (nid)
     1365#ifdef NC_DOUBLE
     1366  ierr = NF_DEF_VAR (nid, "longitude", NF_DOUBLE, 1, idim2,nvarid)
     1367#else
     1368  ierr = NF_DEF_VAR (nid, "longitude", NF_FLOAT, 1, idim2,nvarid)
     1369#endif
     1370  ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 32,  &
     1371     &               "Longitudes de la grille physique")
     1372  ierr = NF_ENDDEF(nid)
     1373#ifdef NC_DOUBLE
     1374  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlon)
     1375#else
     1376  ierr = NF_PUT_VAR_REAL (nid,nvarid,rlon)
     1377#endif
     1378!
     1379  ierr = NF_REDEF (nid)
     1380#ifdef NC_DOUBLE
     1381  ierr = NF_DEF_VAR (nid, "latitude", NF_DOUBLE, 1, idim2,nvarid)
     1382#else
     1383  ierr = NF_DEF_VAR (nid, "latitude", NF_FLOAT, 1, idim2,nvarid)
     1384#endif
     1385  ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 31,    &
     1386     &                        "Latitudes de la grille physique")
     1387  ierr = NF_ENDDEF(nid)
     1388#ifdef NC_DOUBLE
     1389  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlat)
     1390#else
     1391  ierr = NF_PUT_VAR_REAL (nid,nvarid,rlat)
     1392#endif
     1393        ierr = NF_REDEF (nid)
     1394#ifdef NC_DOUBLE
     1395        ierr = NF_DEF_VAR (nid, "TS", NF_DOUBLE, 1, idim2,nvarid)
     1396#else
     1397        ierr = NF_DEF_VAR (nid, "TS", NF_FLOAT, 1, idim2,nvarid)
     1398#endif
     1399        ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 22, &
     1400     &                        "Temperature de surface")
     1401        ierr = NF_ENDDEF(nid)
     1402
     1403
     1404
     1405
     1406  END SUBROUTINE sol_dem_write
     1407!
     1408!#########################################################################
     1409!
     1410  SUBROUTINE gath2cpl(champ_in, champ_out, klon, knon, iim, jjm, knindex)
     1411
     1412! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
     1413! au coupleur.
     1414!
     1415!
     1416! input:         
     1417!   champ_in     champ sur la grille gathere       
     1418!   knon         nombre de points dans le domaine a traiter
     1419!   knindex      index des points de la surface a traiter
     1420!   klon         taille de la grille
     1421!   iim,jjm      dimension de la grille 2D
     1422!
     1423! output:
     1424!   champ_out    champ sur la grille 2D
     1425!
     1426! input
     1427  integer                   :: klon, knon, iim, jjm
     1428  real, dimension(knon)     :: champ_in
     1429  integer, dimension(knon)  :: knindex
     1430! output
     1431  real, dimension(iim,jjm+1)  :: champ_out
     1432! local
     1433  integer                   :: i, ig, j
     1434  real, dimension(klon)     :: tamp
     1435
     1436  do i = 1, knon
     1437    ig = knindex(i)
     1438    tamp(ig) = champ_in(i)
     1439  enddo   
     1440  champ_out(:,1) = tamp(1)
     1441  do j = 2, jjm
     1442    do i = 1, iim
     1443      champ_out(i,j) = tamp((j-2)*jjm + i + 1)
     1444    enddo
     1445  enddo
     1446  champ_out(:,jjm+1) = tamp(klon)
     1447
     1448  END SUBROUTINE gath2cpl
     1449!
     1450!#########################################################################
     1451!
     1452  SUBROUTINE cpl2gath(champ_in, champ_out, klon, knon, iim, jjm, knindex)
     1453
     1454! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
     1455! au coupleur.
     1456!
     1457!
     1458! input:         
     1459!   champ_in     champ sur la grille gathere       
     1460!   knon         nombre de points dans le domaine a traiter
     1461!   knindex      index des points de la surface a traiter
     1462!   klon         taille de la grille
     1463!   iim,jjm      dimension de la grille 2D
     1464!
     1465! output:
     1466!   champ_out    champ sur la grille 2D
     1467!
     1468! input
     1469  integer                   :: klon, knon, iim, jjm
     1470  real, dimension(iim,jjm+1)     :: champ_in
     1471  integer, dimension(knon)  :: knindex
     1472! output
     1473  real, dimension(knon)  :: champ_out
     1474! local
     1475  integer                   :: i, ig, j
     1476  real, dimension(klon)     :: tamp
     1477
     1478  tamp(1) = champ_in(1,1)
     1479  do j = 2, jjm
     1480    do i = 1, iim
     1481      tamp((j-2)*jjm + i + 1) = champ_in(i,j)
     1482    enddo
     1483  enddo
     1484  tamp(klon) = champ_in(1,jjm+1)
     1485
     1486  do i = 1, knon
     1487    ig = knindex(i)
     1488    champ_out(i) = tamp(ig)
     1489  enddo   
     1490
     1491  END SUBROUTINE cpl2gath
     1492!
     1493!#########################################################################
     1494!
    9801495  END MODULE interface_surf
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.F

    r79 r98  
    386386      END
    387387
    388       SUBROUTINE fromcpl(kt, imjm, sst, gla)
     388      SUBROUTINE fromcpl(kt, imjm, sst,sic, alb_sst, alb_sic )
    389389      IMPLICIT none
    390390c
     
    399399      INTEGER imjm, kt
    400400      REAL sst(imjm)          ! sea-surface-temperature
    401       REAL gla(imjm)          ! sea-ice
     401      REAL alb_sst(imjm)  ! open sea albedo
     402      REAL sic(imjm)      ! sea ice cover
     403      REAL alb_sic(imjm)  ! sea ice albedo
     404
    402405c
    403406      INTEGER nuout             ! listing output unit
     
    447450     $          nuout)
    448451            IF (jf.eq.2)
    449      $          CALL locread(cl_read(jf), gla, imjm, nuread, iflag,
     452     $          CALL locread(cl_read(jf), sic, imjm, nuread, iflag,
     453     $          nuout)
     454            IF (jf.eq.3)
     455     $          CALL locread(cl_read(jf), alb_sst, imjm, nuread, iflag,
     456     $          nuout)
     457            IF (jf.eq.4)
     458     $          CALL locread(cl_read(jf), alb_sic, imjm, nuread, iflag,
    450459     $          nuout)
    451460            CLOSE (nuread)
     
    474483c
    475484          CALL SIPC_Read_Model(index, imjm, cmodinf,
    476      $              cljobnam_r,infos, gla)
     485     $              cljobnam_r,infos, sic)
     486c         Index of open sea albedo in total number of fields jpfldo2a:
     487          index = 3
     488c
     489          CALL SIPC_Read_Model(index, imjm, cmodinf,
     490     $              cljobnam_r,infos, alb_sst)
     491c         Index of sea-ice albedo in total number of fields jpfldo2a:
     492          index = 4
     493c
     494          CALL SIPC_Read_Model(index, imjm, cmodinf,
     495     $              cljobnam_r,infos, alb_sic)
    477496c
    478497c
     
    482501c exchanges from ocean=CPL to atmosphere
    483502c
    484           DO jf=1,jpfldo2a
    485             IF (jf.eq.1) CALL CLIM_Import (cl_read(jf) , kt, sst, info)
    486             IF (jf.eq.2) CALL CLIM_Import (cl_read(jf) , kt, gla, info)
    487             IF ( info .NE. CLIM_Ok) THEN
     503        DO jf=1,jpfldo2a
     504          IF (jf.eq.1) CALL CLIM_Import (cl_read(jf) , kt, sst, info)
     505          IF (jf.eq.2) CALL CLIM_Import (cl_read(jf) , kt, sic, info)
     506         IF (jf.eq.3) CALL CLIM_Import (cl_read(jf) , kt, alb_sst, info)
     507         IF (jf.eq.4) CALL CLIM_Import (cl_read(jf) , kt, alb_sic, info)
     508          IF ( info .NE. CLIM_Ok) THEN
    488509                WRITE(nuout,*)'Pb in reading ', cl_read(jf), jf
    489510                WRITE(nuout,*)'Couplage kt is = ',kt
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.h

    r79 r98  
    1       LOGICAL ok_oasis
    2       PARAMETER (ok_oasis = .FALSE.)
    3 c
    4       CHARACTER*4 cchan
    5       PARAMETER (cchan="PIPE")
    6 c      PARAMETER (cchan="CLIM")
     1C
     2C -- oasis.h   
     3C    ******
     4C@
     5C@  Contents : choice for the OASIS version: clim or pipe
     6C@  --------
    77
    8       INTEGER jpmaxfld
    9       PARAMETER(jpmaxfld = 20)
     8      logical ok_oasis
     9      parameter(ok_oasis = .false.)
     10
     11      CHARACTER*8 cchan
     12      PARAMETER ( cchan='PIPE' )
     13C
     14C     --- end of oasis.h
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/param_cou.h

    r79 r98  
    22C -- param_cou.h
    33C
     4        INTEGER jpmaxfld
     5        PARAMETER(jpmaxfld = 100)        ! Number of maximum fields
     6                                         ! exchange betwwen ocean and atmosphere
    47        INTEGER jpflda2o
    5         PARAMETER(jpflda2o = 1)          ! Number of fields exchanged from
     8        PARAMETER(jpflda2o = 8)          ! Number of fields exchanged from
    69                                         ! atmosphere to ocean
    710C
    811        INTEGER jpfldo2a
    9         PARAMETER(jpfldo2a = 1)          ! Number of fields exchanged from
     12        PARAMETER(jpfldo2a = 2)          ! Number of fields exchanged from
    1013                                         ! ocean to atmosphere
    1114C
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/param_sipc.h

    r79 r98  
    1010        INTEGER jpbyteint,jpbyterea, jpbytecha 
    1111        PARAMETER (jpbyteint = 4)
    12         PARAMETER (jpbyterea = 4)
     12        PARAMETER (jpbyterea = 8)
    1313        PARAMETER (jpbytecha = 1)
    1414C@
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/phyetat0.F

    r79 r98  
    11      SUBROUTINE phyetat0 (fichnom,dtime,co2_ppm,solaire,
    2      .            rlat,rlon,tsol,tsoil,deltat,qsol,snow,
     2     .            rlat,rlon, pctsrf, tsol,tsoil,deltat,qsol,snow,
     3     .           albe, evap, rain_fall, snow_fall, solsw, sollw,
    34     .           radsol,rugmer,agesno,clesphy0,
    45     .           zmea,zstd,zsig,zgam,zthe,zpic,zval,rugsrel,tabcntr0,
     
    2728      REAL qsol(klon,nbsrf)
    2829      REAL snow(klon,nbsrf)
     30      REAL albe(klon,nbsrf)
     31      REAL evap(klon,nbsrf)
    2932      REAL radsol(klon)
     33      REAL rain_fall(klon)
     34      REAL snow_fall(klon)
     35      REAL sollw(klon)
     36      real solsw(klon)
    3037      REAL rugmer(klon)
    3138      REAL agesno(klon)
     
    3845      REAL zval(klon)
    3946      REAL rugsrel(klon)
     47      REAL pctsrf(klon, nbsrf)
     48      REAL fractint(klon)
    4049
    4150      REAL t_ancien(klon,klev), q_ancien(klon,klev)
     
    5867c Ouvrir le fichier contenant l'etat initial:
    5968c
     69      print*,'fichnom',fichnom
    6070      ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
    6171      IF (ierr.NE.NF_NOERR) THEN
     
    184194         CALL abort
    185195      ENDIF
    186 c
     196C
     197C
     198C Lecture du masque terre mer
     199C
     200      ierr = NF_INQ_VARID (nid, "masque", nvarid)
     201      IF (ierr .EQ.  NF_NOERR) THEN
     202#ifdef NC_DOUBLE
     203          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zmasq)
     204#else
     205          ierr = NF_GET_VAR_REAL(nid, nvarid, zmasq)
     206#endif
     207          IF (ierr.NE.NF_NOERR) THEN
     208              PRINT*, 'phyetat0: Lecture echouee pour <masque>'
     209              CALL abort
     210          ENDIF
     211      else
     212          PRINT*, 'phyetat0: Le champ <masque> est absent'
     213          PRINT*, 'fichier startphy non compatible avec phyetat0'
     214C      CALL abort
     215      ENDIF
     216C Lecture des fractions pour chaque sous-surface
     217C
     218C initialisation des sous-surfaces
     219C
     220      pctsrf = 0.
     221C
     222C fraction de terre
     223C
     224      ierr = NF_INQ_VARID (nid, "FTER", nvarid)
     225      IF (ierr .EQ.  NF_NOERR) THEN
     226#ifdef NC_DOUBLE
     227          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, pctsrf(1 : klon,is_ter))
     228#else
     229          ierr = NF_GET_VAR_REAL(nid, nvarid, pctsrf(1 : klon,is_ter))
     230#endif
     231          IF (ierr.NE.NF_NOERR) THEN
     232              PRINT*, 'phyetat0: Lecture echouee pour <FTER>'
     233              CALL abort
     234          ENDIF
     235      else
     236          PRINT*, 'phyetat0: Le champ <FTER> est absent'
     237c$$$         CALL abort
     238      ENDIF
     239C
     240C fraction de glace de terre
     241C
     242      ierr = NF_INQ_VARID (nid, "FLIC", nvarid)
     243      IF (ierr .EQ.  NF_NOERR) THEN
     244#ifdef NC_DOUBLE
     245          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, pctsrf(1 : klon,is_lic))
     246#else
     247          ierr = NF_GET_VAR_REAL(nid, nvarid, pctsrf(1 : klon,is_lic))
     248#endif
     249          IF (ierr.NE.NF_NOERR) THEN
     250              PRINT*, 'phyetat0: Lecture echouee pour <FLIC>'
     251              CALL abort
     252          ENDIF
     253      else
     254          PRINT*, 'phyetat0: Le champ <FLIC> est absent'
     255c$$$         CALL abort
     256      ENDIF
     257C
     258C fraction d'ocean
     259C
     260      ierr = NF_INQ_VARID (nid, "FOCE", nvarid)
     261      IF (ierr .EQ.  NF_NOERR) THEN
     262#ifdef NC_DOUBLE
     263          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, pctsrf(1 : klon,is_oce))
     264#else
     265          ierr = NF_GET_VAR_REAL(nid, nvarid, pctsrf(1 : klon,is_oce))
     266#endif
     267          IF (ierr.NE.NF_NOERR) THEN
     268              PRINT*, 'phyetat0: Lecture echouee pour <FOCE>'
     269              CALL abort
     270          ENDIF
     271      else
     272          PRINT*, 'phyetat0: Le champ <FOCE> est absent'
     273c$$$         CALL abort
     274      ENDIF
     275C
     276C fraction glace de mer
     277C
     278      ierr = NF_INQ_VARID (nid, "FSIC", nvarid)
     279      IF (ierr .EQ.  NF_NOERR) THEN
     280#ifdef NC_DOUBLE
     281          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, pctsrf(1 : klon,is_sic))
     282#else
     283          ierr = NF_GET_VAR_REAL(nid, nvarid, pctsrf(1 : klon, is_sic))
     284#endif
     285          IF (ierr.NE.NF_NOERR) THEN
     286              PRINT*, 'phyetat0: Lecture echouee pour <FSIC>'
     287              CALL abort
     288          ENDIF
     289      else
     290          PRINT*, 'phyetat0: Le champ <FSIC> est absent'
     291c$$$         CALL abort
     292      ENDIF
     293C
     294C  Verification de l'adequation entre le masque et les sous-surfaces
     295C
     296      fractint( 1 : klon) = pctsrf(1 : klon, is_ter)
     297     $    + pctsrf(1 : klon, is_lic)
     298      DO i = 1 , klon
     299        IF ( abs(fractint(i) - zmasq(i) ) .GT. EPSFRA ) THEN
     300            WRITE(*,*) 'phyetat0: attention fraction terre pas ',
     301     $          'coherente ', i, zmasq(i), pctsrf(i, is_ter)
     302     $          ,pctsrf(i, is_lic)
     303        ENDIF
     304      END DO
     305      fractint (1 : klon) =  pctsrf(1 : klon, is_oce)
     306     $    + pctsrf(1 : klon, is_sic)
     307      DO i = 1 , klon
     308        IF ( abs( fractint(i) - (1. - zmasq(i))) .GT. EPSFRA ) THEN
     309            WRITE(*,*) 'phyetat0 attention fraction ocean pas ',
     310     $          'coherente ', i, zmasq(i) , pctsrf(i, is_oce)
     311     $          ,pctsrf(i, is_sic)
     312        ENDIF
     313      END DO
     314C
    187315c Lecture des temperatures du sol:
    188316c
     
    419547      ENDIF
    420548c
     549c Lecture de albedo au sol:
     550c
     551      ierr = NF_INQ_VARID (nid, "ALBE", nvarid)
     552      IF (ierr.NE.NF_NOERR) THEN
     553         PRINT*, 'phyetat0: Le champ <ALBE> est absent'
     554         PRINT*, '          Mais je vais essayer de lire ALBE**'
     555         DO nsrf = 1, nbsrf
     556           IF (nsrf.GT.99) THEN
     557             PRINT*, "Trop de sous-mailles"
     558             CALL abort
     559           ENDIF
     560           WRITE(str2,'(i2.2)') nsrf
     561           ierr = NF_INQ_VARID (nid, "ALBE"//str2, nvarid)
     562           IF (ierr.NE.NF_NOERR) THEN
     563              PRINT*, "phyetat0: Le champ <ALBE"//str2//"> est absent"
     564              CALL abort
     565           ENDIF
     566#ifdef NC_DOUBLE
     567           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, albe(1,nsrf))
     568#else
     569           ierr = NF_GET_VAR_REAL(nid, nvarid, albe(1,nsrf))
     570#endif
     571           IF (ierr.NE.NF_NOERR) THEN
     572             PRINT*, "phyetat0: Lecture echouee pour <ALBE"//str2//">"
     573             CALL abort
     574           ENDIF
     575           xmin = 1.0E+20
     576           xmax = -1.0E+20
     577           DO i = 1, klon
     578              xmin = MIN(snow(i,nsrf),xmin)
     579              xmax = MAX(snow(i,nsrf),xmax)
     580           ENDDO
     581           PRINT*,'Neige du sol ALBE**:', nsrf, xmin, xmax
     582         ENDDO
     583      ELSE
     584         PRINT*, 'phyetat0: Le champ <ALBE> est present'
     585         PRINT*, '          J ignore donc les autres ALBE**'
     586#ifdef NC_DOUBLE
     587         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, albe(1,1))
     588#else
     589         ierr = NF_GET_VAR_REAL(nid, nvarid, albe(1,1))
     590#endif
     591         IF (ierr.NE.NF_NOERR) THEN
     592            PRINT*, "phyetat0: Lecture echouee pour <ALBE>"
     593            CALL abort
     594         ENDIF
     595         xmin = 1.0E+20
     596         xmax = -1.0E+20
     597         DO i = 1, klon
     598            xmin = MIN(albe(i,1),xmin)
     599            xmax = MAX(albe(i,1),xmax)
     600         ENDDO
     601         PRINT*,'Neige du sol <ALBE>', xmin, xmax
     602         DO nsrf = 2, nbsrf
     603         DO i = 1, klon
     604            albe(i,nsrf) = albe(i,1)
     605         ENDDO
     606         ENDDO
     607      ENDIF
     608
     609c
     610c Lecture de evaporation: 
     611c
     612      ierr = NF_INQ_VARID (nid, "EVAP", nvarid)
     613      IF (ierr.NE.NF_NOERR) THEN
     614         PRINT*, 'phyetat0: Le champ <EVAP> est absent'
     615         PRINT*, '          Mais je vais essayer de lire EVAP**'
     616         DO nsrf = 1, nbsrf
     617           IF (nsrf.GT.99) THEN
     618             PRINT*, "Trop de sous-mailles"
     619             CALL abort
     620           ENDIF
     621           WRITE(str2,'(i2.2)') nsrf
     622           ierr = NF_INQ_VARID (nid, "EVAP"//str2, nvarid)
     623           IF (ierr.NE.NF_NOERR) THEN
     624              PRINT*, "phyetat0: Le champ <EVAP"//str2//"> est absent"
     625              CALL abort
     626           ENDIF
     627#ifdef NC_DOUBLE
     628           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, evap(1,nsrf))
     629#else
     630           ierr = NF_GET_VAR_REAL(nid, nvarid, evap(1,nsrf))
     631#endif
     632           IF (ierr.NE.NF_NOERR) THEN
     633             PRINT*, "phyetat0: Lecture echouee pour <EVAP"//str2//">"
     634             CALL abort
     635           ENDIF
     636           xmin = 1.0E+20
     637           xmax = -1.0E+20
     638           DO i = 1, klon
     639              xmin = MIN(evap(i,nsrf),xmin)
     640              xmax = MAX(evap(i,nsrf),xmax)
     641           ENDDO
     642           PRINT*,'Neige du sol EVAP**:', nsrf, xmin, xmax
     643         ENDDO
     644      ELSE
     645         PRINT*, 'phyetat0: Le champ <EVAP> est present'
     646         PRINT*, '          J ignore donc les autres EVAP**'
     647#ifdef NC_DOUBLE
     648         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, evap(1,1))
     649#else
     650         ierr = NF_GET_VAR_REAL(nid, nvarid, evap(1,1))
     651#endif
     652         IF (ierr.NE.NF_NOERR) THEN
     653            PRINT*, "phyetat0: Lecture echouee pour <EVAP>"
     654            CALL abort
     655         ENDIF
     656         xmin = 1.0E+20
     657         xmax = -1.0E+20
     658         DO i = 1, klon
     659            xmin = MIN(evap(i,1),xmin)
     660            xmax = MAX(evap(i,1),xmax)
     661         ENDDO
     662         PRINT*,'Neige du sol <EVAP>', xmin, xmax
     663         DO nsrf = 2, nbsrf
     664         DO i = 1, klon
     665            evap(i,nsrf) = evap(i,1)
     666         ENDDO
     667         ENDDO
     668      ENDIF
     669c
     670c Lecture precipitation liquide:
     671c
     672      ierr = NF_INQ_VARID (nid, "rain_f", nvarid)
     673      IF (ierr.NE.NF_NOERR) THEN
     674         PRINT*, 'phyetat0: Le champ <rain_f> est absent'
     675         CALL abort
     676      ENDIF
     677#ifdef NC_DOUBLE
     678      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rain_fall)
     679#else
     680      ierr = NF_GET_VAR_REAL(nid, nvarid, rain_fall)
     681#endif
     682      IF (ierr.NE.NF_NOERR) THEN
     683         PRINT*, 'phyetat0: Lecture echouee pour <rain_f>'
     684         CALL abort
     685      ENDIF
     686      xmin = 1.0E+20
     687      xmax = -1.0E+20
     688      DO i = 1, klon
     689         xmin = MIN(rain_fall(i),xmin)
     690         xmax = MAX(rain_fall(i),xmax)
     691      ENDDO
     692      PRINT*,'Precipitation liquide rain_f:', xmin, xmax
     693c
     694c Lecture precipitation solide:
     695c
     696      ierr = NF_INQ_VARID (nid, "snow_f", nvarid)
     697      IF (ierr.NE.NF_NOERR) THEN
     698         PRINT*, 'phyetat0: Le champ <snow_f> est absent'
     699         CALL abort
     700      ENDIF
     701#ifdef NC_DOUBLE
     702      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, snow_fall)
     703#else
     704      ierr = NF_GET_VAR_REAL(nid, nvarid, snow_fall)
     705#endif
     706      IF (ierr.NE.NF_NOERR) THEN
     707         PRINT*, 'phyetat0: Lecture echouee pour <snow_f>'
     708         CALL abort
     709      ENDIF
     710      xmin = 1.0E+20
     711      xmax = -1.0E+20
     712      DO i = 1, klon
     713         xmin = MIN(snow_fall(i),xmin)
     714         xmax = MAX(snow_fall(i),xmax)
     715      ENDDO
     716      PRINT*,'Precipitation solide snow_f:', xmin, xmax
     717c
     718c Lecture rayonnement solaire au sol:
     719c
     720      ierr = NF_INQ_VARID (nid, "solsw", nvarid)
     721      IF (ierr.NE.NF_NOERR) THEN
     722         PRINT*, 'phyetat0: Le champ <solsw> est absent'
     723         PRINT*, 'mis a zero'
     724         solsw = 0.
     725      ELSE
     726#ifdef NC_DOUBLE
     727        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, solsw)
     728#else
     729        ierr = NF_GET_VAR_REAL(nid, nvarid, solsw)
     730#endif
     731        IF (ierr.NE.NF_NOERR) THEN
     732          PRINT*, 'phyetat0: Lecture echouee pour <solsw>'
     733          CALL abort
     734        ENDIF
     735      ENDIF
     736      xmin = 1.0E+20
     737      xmax = -1.0E+20
     738      DO i = 1, klon
     739         xmin = MIN(solsw(i),xmin)
     740         xmax = MAX(solsw(i),xmax)
     741      ENDDO
     742      PRINT*,'Rayonnement solaire au sol solsw:', xmin, xmax
     743c
     744c Lecture rayonnement IF au sol:
     745c
     746      ierr = NF_INQ_VARID (nid, "sollw", nvarid)
     747      IF (ierr.NE.NF_NOERR) THEN
     748         PRINT*, 'phyetat0: Le champ <sollw> est absent'
     749         PRINT*, 'mis a zero'
     750         sollw = 0.
     751      ELSE
     752#ifdef NC_DOUBLE
     753        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, sollw)
     754#else
     755        ierr = NF_GET_VAR_REAL(nid, nvarid, sollw)
     756#endif
     757        IF (ierr.NE.NF_NOERR) THEN
     758          PRINT*, 'phyetat0: Lecture echouee pour <sollw>'
     759          CALL abort
     760        ENDIF
     761      ENDIF
     762      xmin = 1.0E+20
     763      xmax = -1.0E+20
     764      DO i = 1, klon
     765         xmin = MIN(sollw(i),xmin)
     766         xmax = MAX(sollw(i),xmax)
     767      ENDDO
     768      PRINT*,'Rayonnement IF au sol sollw:', xmin, xmax
     769
     770c
    421771c Lecture du rayonnement net au sol:
    422772c
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/phyredem.F

    r79 r98  
    11      SUBROUTINE phyredem (fichnom,dtime,radpas,co2_ppm,solaire,
    2      .           rlat,rlon,tsol,tsoil,deltat,qsol,snow,
     2     .           rlat,rlon, pctsrf,tsol,tsoil,deltat,qsol,snow,
     3     .           albedo, evap, rain_fall, snow_fall,
     4     .           solsw, sollw,
    35     .           radsol,rugmer,agesno,
    46     .           zmea,zstd,zsig,zgam,zthe,zpic,zval,rugsrel,
     
    2931      REAL qsol(klon,nbsrf)
    3032      REAL snow(klon,nbsrf)
     33      REAL albedo(klon,nbsrf)
     34      REAL evap(klon,nbsrf)
     35      REAL rain_fall(klon)
     36      REAL snow_fall(klon)
     37      real solsw(klon)
     38      real sollw(klon)
    3139      REAL radsol(klon)
    3240      REAL rugmer(klon)
     
    4048      REAL zval(klon)
    4149      REAL rugsrel(klon)
     50      REAL pctsrf(klon, nbsrf)
    4251      REAL t_ancien(klon,klev), q_ancien(klon,klev)
    4352c
     
    134143#endif
    135144c
     145C PB ajout du masque terre/mer
     146C
     147      ierr = NF_REDEF (nid)
     148#ifdef NC_DOUBLE
     149      ierr = NF_DEF_VAR (nid, "masque", NF_DOUBLE, 1, idim2,nvarid)
     150#else
     151      ierr = NF_DEF_VAR (nid, "masque", NF_FLOAT, 1, idim2,nvarid)
     152#endif
     153      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 16,
     154     .                        "masque terre mer")
     155      ierr = NF_ENDDEF(nid)
     156#ifdef NC_DOUBLE
     157      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zmasq)
     158#else
     159      ierr = NF_PUT_VAR_REAL (nid,nvarid,zmasq)
     160#endif     
     161c BP ajout des fraction de chaque sous-surface
     162C
     163C 1. fraction de terre
     164C
     165      ierr = NF_REDEF (nid)
     166#ifdef NC_DOUBLE
     167      ierr = NF_DEF_VAR (nid, "FTER", NF_DOUBLE, 1, idim2,nvarid)
     168#else
     169      ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 1, idim2,nvarid)
     170#endif
     171      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 21,
     172     .                        "fraction de continent")
     173      ierr = NF_ENDDEF(nid)
     174#ifdef NC_DOUBLE
     175      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pctsrf(1 : klon, is_ter))
     176#else
     177      ierr = NF_PUT_VAR_REAL (nid,nvarid,pctsrf(1 : klon, is_ter))
     178#endif
     179C
     180C 2. Fraction de glace de terre
     181C
     182      ierr = NF_REDEF (nid)
     183#ifdef NC_DOUBLE
     184      ierr = NF_DEF_VAR (nid, "FLIC", NF_DOUBLE, 1, idim2,nvarid)
     185#else
     186      ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 1, idim2,nvarid)
     187#endif
     188      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 24,
     189     .                        "fraction glace de terre")
     190      ierr = NF_ENDDEF(nid)
     191#ifdef NC_DOUBLE
     192      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pctsrf(1 : klon,is_lic))
     193#else
     194      ierr = NF_PUT_VAR_REAL (nid,nvarid,pctsrf(1 : klon, is_lic))
     195#endif
     196C
     197C 3. fraction ocean
     198C
     199      ierr = NF_REDEF (nid)
     200#ifdef NC_DOUBLE
     201      ierr = NF_DEF_VAR (nid, "FOCE", NF_DOUBLE, 1, idim2,nvarid)
     202#else
     203      ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 1, idim2,nvarid)
     204#endif
     205      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 14,
     206     .                        "fraction ocean")
     207      ierr = NF_ENDDEF(nid)
     208#ifdef NC_DOUBLE
     209      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pctsrf(1 : klon, is_oce))
     210#else
     211      ierr = NF_PUT_VAR_REAL (nid,nvarid,pctsrf(1 : klon, is_oce))
     212#endif
     213C
     214C 4. Fraction glace de mer
     215C
     216      ierr = NF_REDEF (nid)
     217#ifdef NC_DOUBLE
     218      ierr = NF_DEF_VAR (nid, "FSIC", NF_DOUBLE, 1, idim2,nvarid)
     219#else
     220      ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 1, idim2,nvarid)
     221#endif
     222      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 18,
     223     .                        "fraction glace mer")
     224      ierr = NF_ENDDEF(nid)
     225#ifdef NC_DOUBLE
     226      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pctsrf(1 : klon, is_sic))
     227#else
     228      ierr = NF_PUT_VAR_REAL (nid,nvarid,pctsrf(1 : klon, is_sic))
     229#endif
     230C
     231C
    136232c
    137233      DO nsrf = 1, nbsrf
     
    227323        ierr = NF_REDEF (nid)
    228324#ifdef NC_DOUBLE
     325        ierr = NF_DEF_VAR (nid,"ALBE"//str2,NF_DOUBLE,1,idim2,nvarid)
     326#else
     327        ierr = NF_DEF_VAR (nid,"ALBE"//str2,NF_FLOAT,1,idim2,nvarid)
     328#endif
     329        ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 23,
     330     .                        "albedo de surface No."//str2)
     331        ierr = NF_ENDDEF(nid)
     332        ELSE
     333        PRINT*, "Trop de sous-mailles"
     334        CALL abort
     335        ENDIF
     336#ifdef NC_DOUBLE
     337      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,albedo(1,nsrf))
     338#else
     339      ierr = NF_PUT_VAR_REAL (nid,nvarid,albedo(1,nsrf))
     340#endif
     341      ENDDO
     342c
     343      DO nsrf = 1, nbsrf
     344        IF (nsrf.LE.99) THEN
     345        WRITE(str2,'(i2.2)') nsrf
     346        ierr = NF_REDEF (nid)
     347#ifdef NC_DOUBLE
     348        ierr = NF_DEF_VAR (nid,"EVAP"//str2,NF_DOUBLE,1,idim2,nvarid)
     349#else
     350        ierr = NF_DEF_VAR (nid,"EVAP"//str2,NF_FLOAT,1,idim2,nvarid)
     351#endif
     352        ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
     353     .                        "Evaporation de surface No."//str2)
     354        ierr = NF_ENDDEF(nid)
     355        ELSE
     356        PRINT*, "Trop de sous-mailles"
     357        CALL abort
     358        ENDIF
     359#ifdef NC_DOUBLE
     360      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,evap(1,nsrf))
     361#else
     362      ierr = NF_PUT_VAR_REAL (nid,nvarid,evap(1,nsrf))
     363#endif
     364      ENDDO
     365
     366c
     367      DO nsrf = 1, nbsrf
     368        IF (nsrf.LE.99) THEN
     369        WRITE(str2,'(i2.2)') nsrf
     370        ierr = NF_REDEF (nid)
     371#ifdef NC_DOUBLE
    229372        ierr = NF_DEF_VAR (nid,"SNOW"//str2,NF_DOUBLE,1,idim2,nvarid)
    230373#else
     
    244387#endif
    245388      ENDDO
     389
    246390c
    247391      ierr = NF_REDEF (nid)
     
    262406      ierr = NF_REDEF (nid)
    263407#ifdef NC_DOUBLE
     408      ierr = NF_DEF_VAR (nid, "solsw", NF_DOUBLE, 1, idim2,nvarid)
     409#else
     410      ierr = NF_DEF_VAR (nid, "solsw", NF_FLOAT, 1, idim2,nvarid)
     411#endif
     412      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 32,
     413     .                        "Rayonnement solaire a la surface")
     414      ierr = NF_ENDDEF(nid)
     415#ifdef NC_DOUBLE
     416      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,solsw)
     417#else
     418      ierr = NF_PUT_VAR_REAL (nid,nvarid,solsw)
     419#endif
     420c
     421      ierr = NF_REDEF (nid)
     422#ifdef NC_DOUBLE
     423      ierr = NF_DEF_VAR (nid, "sollw", NF_DOUBLE, 1, idim2,nvarid)
     424#else
     425      ierr = NF_DEF_VAR (nid, "sollw", NF_FLOAT, 1, idim2,nvarid)
     426#endif
     427      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 27,
     428     .                        "Rayonnement IF a la surface")
     429      ierr = NF_ENDDEF(nid)
     430#ifdef NC_DOUBLE
     431      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,sollw)
     432#else
     433      ierr = NF_PUT_VAR_REAL (nid,nvarid,sollw)
     434#endif
     435c
     436      ierr = NF_REDEF (nid)
     437#ifdef NC_DOUBLE
     438      ierr = NF_DEF_VAR (nid, "rain_f", NF_DOUBLE, 1, idim2,nvarid)
     439#else
     440      ierr = NF_DEF_VAR (nid, "rain_f", NF_FLOAT, 1, idim2,nvarid)
     441#endif
     442      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 21,
     443     .                        "precipitation liquide")
     444      ierr = NF_ENDDEF(nid)
     445#ifdef NC_DOUBLE
     446      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rain_fall)
     447#else
     448      ierr = NF_PUT_VAR_REAL (nid,nvarid,rain_fall)
     449#endif
     450c
     451      ierr = NF_REDEF (nid)
     452#ifdef NC_DOUBLE
     453      ierr = NF_DEF_VAR (nid, "snow_f", NF_DOUBLE, 1, idim2,nvarid)
     454#else
     455      ierr = NF_DEF_VAR (nid, "snow_f", NF_FLOAT, 1, idim2,nvarid)
     456#endif
     457      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
     458     .                        "precipitation solide")
     459      ierr = NF_ENDDEF(nid)
     460#ifdef NC_DOUBLE
     461      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,snow_fall)
     462#else
     463      ierr = NF_PUT_VAR_REAL (nid,nvarid,snow_fall)
     464#endif
     465c
     466      ierr = NF_REDEF (nid)
     467#ifdef NC_DOUBLE
    264468      ierr = NF_DEF_VAR (nid, "RUGMER", NF_DOUBLE, 1, idim2,nvarid)
    265469#else
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F

    r88 r98  
    100100      REAL soilcap(klon,nbsrf), soilflux(klon,nbsrf)
    101101      SAVE soilcap, soilflux
     102      logical ok_veget
     103      parameter (ok_veget = .false.)
    102104c======================================================================
    103105c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
     
    132134      INTEGER iliq          ! indice de traceurs pour eau liquide
    133135      PARAMETER (iliq=2)
    134 c
     136
    135137      INTEGER nvm           ! nombre de vegetations
    136138      PARAMETER (nvm=8)
    137139      REAL veget(klon,nvm)  ! couverture vegetale
    138140      SAVE veget
     141
     142c
    139143c
    140144c Variables argument:
     
    224228      SAVE ftsoil                 ! temperature dans le sol
    225229c
     230      REAL fevap(klon,nbsrf)
     231      SAVE fevap                 ! evaporation
     232c
    226233      REAL deltat(klon)
    227234      SAVE deltat                 ! ecart avec la SST de reference
     
    232239      REAL fsnow(klon,nbsrf)
    233240      SAVE fsnow                  ! epaisseur neigeuse
     241c
     242      REAL falbe(klon,nbsrf)
     243      SAVE falbe                  ! albedo par type de surface
    234244c
    235245      REAL rugmer(klon)
     
    379389      REAL cldemi(klon,klev)  ! emissivite infrarouge
    380390c
    381       REAL fluxq(klon,klev)   ! flux turbulent d'humidite
    382       REAL fluxt(klon,klev)   ! flux turbulent de chaleur
    383       REAL fluxu(klon,klev)   ! flux turbulent de vitesse u
    384       REAL fluxv(klon,klev)   ! flux turbulent de vitesse v
    385 c
     391C§§§ PB
     392      REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
     393      REAL fluxt(klon,klev, nbsrf)   ! flux turbulent de chaleur
     394      REAL fluxu(klon,klev, nbsrf)   ! flux turbulent de vitesse u
     395      REAL fluxv(klon,klev, nbsrf)   ! flux turbulent de vitesse v
     396c
     397      REAL zxfluxt(klon, klev)
     398      REAL zxfluxq(klon, klev)
     399      REAL zxfluxu(klon, klev)
     400      REAL zxfluxv(klon, klev)
     401C§§§
    386402      REAL heat(klon,klev)    ! chauffage solaire
    387403      REAL heat0(klon,klev)   ! chauffage solaire ciel clair
     
    424440c
    425441      REAL zphi(klon,klev)
    426       REAL zx_tmp_x(iim), zx_tmp_yjjmp1
    427442      REAL zx_relief(iim,jjmp1)
    428443      REAL zx_aire(iim,jjmp1)
     
    561576c
    562577       IF (debut) THEN
    563 c
    564  
    565          IF (ok_oasis) THEN
    566             PRINT*, "Attentions! les parametres suivants sont fixes:"
    567             PRINT *,'***********************************************'
    568             PRINT*, "npas, nexca, itimestep=", npas, nexca, itimestep
    569             PRINT*, "Changer-les manuellement s il le faut"
    570             PRINT *,'***********************************************'
    571             CALL inicma( npas, nexca, itimestep)
    572          ENDIF
    573 c
    574          IF (ok_ocean) THEN
    575             PRINT*, '************************'
    576             PRINT*, 'SLAB OCEAN est active, prenez precautions !'
    577             PRINT*, '************************'
    578          ENDIF
    579 c
     578
    580579         DO k = 2, nvm          ! pas de vegetation
    581580            DO i = 1, klon
     
    588587         PRINT*, 'Pas de vegetation; desert partout'
    589588c
     589c
    590590c Initialiser les compteurs:
    591591c
     
    595595c
    596596         CALL phyetat0 ("startphy.nc",dtime,co2_ppm,solaire,
    597      .         rlat,rlon,ftsol,ftsoil,deltat,fqsol,fsnow,
    598      .        radsol,rugmer,agesno,clesphy0,
     597     .       rlat,rlon,pctsrf, ftsol,ftsoil,deltat,fqsol,fsnow,
     598     .       falbe, fevap, rain_fall,snow_fall,sollw, solsw,
     599     .       radsol,rugmer,agesno,clesphy0,
    599600     .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
    600601     .       t_ancien, q_ancien, ancien_ok )
     
    646647         ENDIF
    647648c
    648          IF (soil_model) THEN
    649             DO nsrf = 1, nbsrf
    650             CALL soil(dtime, nsrf, fsnow(1,nsrf),
    651      .                ftsol(1,nsrf), ftsoil(1,1,nsrf),
    652      .                soilcap(1,nsrf), soilflux(1,nsrf))
    653             ENDDO
    654          ENDIF
    655649c
    656650         lmt_pas = NINT(86400./dtime * 1.0)   ! tous les jours
     
    716710     .                "ave(X)", zsto,zout)
    717711c
     712         CALL histdef(nid_day, "tter", "Surface Temperature", "K",
     713     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     714     .                "ave(X)", zsto,zout)
     715c
     716         CALL histdef(nid_day, "tlic", "Surface Temperature", "K",
     717     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     718     .                "ave(X)", zsto,zout)
     719c
     720         CALL histdef(nid_day, "toce", "Surface Temperature", "K",
     721     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     722     .                "ave(X)", zsto,zout)
     723c
     724         CALL histdef(nid_day, "tsic", "Surface Temperature", "K",
     725     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     726     .                "ave(X)", zsto,zout)
     727c
    718728         CALL histdef(nid_day, "psol", "Surface Pressure", "Pa",
    719729     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     
    768778     .                "ave(X)", zsto,zout)
    769779c
     780C §§§ PB flux pour chauqe sous surface
     781C
     782         DO nsrf = 1, nbsrf
     783C
     784           call histdef(nid_day, "pourc_"//clnsurf(nsrf),
     785     $         "Fraction"//clnsurf(nsrf), "W/m2", 
     786     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
     787     $         "ave(X)", zsto,zout)
     788
     789           call histdef(nid_day, "sens_"//clnsurf(nsrf),
     790     $         "Sensible heat flux "//clnsurf(nsrf), "W/m2", 
     791     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
     792     $         "ave(X)", zsto,zout)
     793c
     794           call histdef(nid_day, "lat_"//clnsurf(nsrf),
     795     $         "Latent heat flux "//clnsurf(nsrf), "W/m2", 
     796     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
     797     $         "ave(X)", zsto,zout)
     798C
     799           call histdef(nid_day, "taux_"//clnsurf(nsrf),
     800     $         "Zonal wind stress"//clnsurf(nsrf),"Pa",
     801     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
     802     $         "ave(X)", zsto,zout)
     803
     804           call histdef(nid_day, "tauy_"//clnsurf(nsrf),
     805     $         "Meridional xind stress "//clnsurf(nsrf), "Pa", 
     806     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
     807     $         "ave(X)", zsto,zout)
     808C§§§
     809         END DO
     810           
    770811         CALL histdef(nid_day, "ruis", "Runoff", "mm/day",
    771812     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     
    956997     .                "ave(X)", zsto,zout)
    957998c
     999         DO nsrf = 1, nbsrf
     1000C
     1001           call histdef(nid_mth, "pourc_"//clnsurf(nsrf),
     1002     $         "Fraction "//clnsurf(nsrf), "W/m2", 
     1003     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
     1004     $         "ave(X)", zsto,zout)
     1005C
     1006           call histdef(nid_mth, "sens_"//clnsurf(nsrf),
     1007     $         "Sensible heat flux "//clnsurf(nsrf), "W/m2", 
     1008     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
     1009     $         "ave(X)", zsto,zout)
     1010c
     1011           call histdef(nid_mth, "lat_"//clnsurf(nsrf),
     1012     $         "Latent heat flux "//clnsurf(nsrf), "W/m2", 
     1013     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
     1014     $         "ave(X)", zsto,zout)
     1015C
     1016           call histdef(nid_mth, "taux_"//clnsurf(nsrf),
     1017     $         "Zonal wind stress"//clnsurf(nsrf), "Pa", 
     1018     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
     1019     $         "ave(X)", zsto,zout)
     1020
     1021           call histdef(nid_mth, "tauy_"//clnsurf(nsrf),
     1022     $         "Meridional xind stress "//clnsurf(nsrf), "Pa", 
     1023     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
     1024     $         "ave(X)", zsto,zout)
     1025         END DO
     1026C
    9581027         CALL histdef(nid_mth, "ruis", "Runoff", "mm/day",
    9591028     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     
    12091278c Champs 2D:
    12101279c
    1211          CALL histdef(nid_ins, "psol", "Surface Pressure", "Pa",
     1280        CALL histdef(nid_ins, "tsol", "Surface Temperature", "K",
     1281     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     1282     .                "inst(X)", zsto,zout)
     1283c
     1284        CALL histdef(nid_ins, "psol", "Surface Pressure", "Pa",
    12121285     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
    12131286     .                "inst(X)", zsto,zout)
     
    12561329     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
    12571330     .                "inst(X)", zsto,zout)
     1331
     1332         DO nsrf = 1, nbsrf
     1333C
     1334           call histdef(nid_ins, "pourc_"//clnsurf(nsrf),
     1335     $         "Fraction"//clnsurf(nsrf), "W/m2", 
     1336     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
     1337     $         "inst(X)", zsto,zout)
     1338
     1339           call histdef(nid_ins, "sens_"//clnsurf(nsrf),
     1340     $         "Sensible heat flux "//clnsurf(nsrf), "W/m2", 
     1341     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
     1342     $         "inst(X)", zsto,zout)
     1343c
     1344           call histdef(nid_ins, "tsol_"//clnsurf(nsrf),
     1345     $         "Surface Temperature"//clnsurf(nsrf), "W/m2", 
     1346     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
     1347     $         "inst(X)", zsto,zout)
     1348c
     1349           call histdef(nid_ins, "lat_"//clnsurf(nsrf),
     1350     $         "Latent heat flux "//clnsurf(nsrf), "W/m2", 
     1351     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
     1352     $         "inst(X)", zsto,zout)
     1353C
     1354           call histdef(nid_ins, "taux_"//clnsurf(nsrf),
     1355     $         "Zonal wind stress"//clnsurf(nsrf),"Pa",
     1356     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
     1357     $         "inst(X)", zsto,zout)
     1358
     1359           call histdef(nid_ins, "tauy_"//clnsurf(nsrf),
     1360     $         "Meridional xind stress "//clnsurf(nsrf), "Pa", 
     1361     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
     1362     $         "inst(X)", zsto,zout)
     1363C§§§
     1364         END DO
    12581365c
    12591366c Champs 3D:
     
    13061413cc         ENDDO
    13071414c
    1308          IF (ok_oasis) THEN
    1309          DO i = 1, klon
    1310            oas_sols(i) = 0.0
    1311            oas_nsol(i) = 0.0
    1312            oas_rain(i) = 0.0
    1313            oas_snow(i) = 0.0
    1314            oas_evap(i) = 0.0
    1315            oas_ruis(i) = 0.0
    1316            oas_tsol(i) = 0.0
    1317            oas_fder(i) = 0.0
    1318            oas_albe(i) = 0.0
    1319            oas_taux(i) = 0.0
    1320            oas_tauy(i) = 0.0
    1321          ENDDO
    1322          ENDIF
    13231415c
    13241416      ENDIF
     
    14211513         CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
    14221514      ENDIF
    1423 cccccccccc
    1424       IF (ok_oasis .AND. MOD(itap-1,nexca).EQ.0) THEN
    1425 C
    1426          CALL fromcpl(itap,jjmp1*iim,
    1427      .        cpl_sst,cpl_sic,cpl_alb_sst,cpl_alb_sic)
    1428          DO i = 1, iim-1 ! un seul point pour le pole nord
    1429             cpl_sst(i,1) = cpl_sst(iim,1)
    1430             cpl_sic(i,1) = cpl_sic(iim,1)
    1431             cpl_alb_sst(i,1) = cpl_alb_sst(iim,1)
    1432             cpl_alb_sic(i,1) = cpl_alb_sic(iim,1)
    1433          ENDDO
    1434          DO i = 2, iim ! un seul point pour le pole sud
    1435             cpl_sst(i,jjmp1) = cpl_sst(1,jjmp1)
    1436             cpl_sic(i,jjmp1) = cpl_sic(1,jjmp1)
    1437             cpl_alb_sst(i,jjmp1) = cpl_alb_sst(1,jjmp1)
    1438             cpl_alb_sic(i,jjmp1) = cpl_alb_sic(1,jjmp1)
    1439          ENDDO
    1440 c
    1441          ig = 1
    1442          IF (pctsrf(ig,is_oce).GT.epsfra .OR.
    1443      .       pctsrf(ig,is_sic).GT.epsfra) THEN
    1444             pctsrf(ig,is_oce) = pctsrf(ig,is_oce)
    1445      .                        - (cpl_sic(1,1)-pctsrf(ig,is_sic))
    1446             pctsrf(ig,is_sic) = cpl_sic(1,1)
    1447             lmt_sst(ig) = cpl_sst(1,1)
    1448          ENDIF
    1449          DO j = 2, jjm
    1450          DO i = 1, iim
    1451          ig = ig + 1
    1452          IF (pctsrf(ig,is_oce).GT.epsfra .OR.
    1453      .       pctsrf(ig,is_sic).GT.epsfra) THEN
    1454            pctsrf(ig,is_oce) = pctsrf(ig,is_oce)
    1455      .                       - (cpl_sic(i,j)-pctsrf(ig,is_sic))
    1456            pctsrf(ig,is_sic) = cpl_sic(i,j)
    1457            lmt_sst(ig) = cpl_sst(i,j)
    1458          ENDIF
    1459          ENDDO
    1460          ENDDO
    1461          ig = ig + 1
    1462          IF (pctsrf(ig,is_oce).GT.epsfra .OR.
    1463      .       pctsrf(ig,is_sic).GT.epsfra) THEN
    1464             pctsrf(ig,is_oce) = pctsrf(ig,is_oce)
    1465      .                        - (cpl_sic(1,jjmp1)-pctsrf(ig,is_sic))
    1466             pctsrf(ig,is_sic) = cpl_sic(1,jjmp1)
    1467             lmt_sst(ig) = cpl_sst(1,jjmp1)
    1468          ENDIF
    1469 c
    1470       ENDIF ! ok_oasis
    1471 cccccccccc
    1472 c
    1473  
    1474       IF (ok_ocean) THEN
    1475          DO i = 1, klon
    1476             ftsol(i,is_oce) = lmt_sst(i) + deltat(i)
    1477          ENDDO
    1478 
    1479       ELSE
    1480          DO i = 1, klon
    1481             ftsol(i,is_oce) = lmt_sst(i)
    1482          ENDDO
    1483 
    1484       ENDIF
    14851515c
    14861516c Re-evaporer l'eau liquide nuageuse
     
    15231553c
    15241554      CALL clmain(dtime,pctsrf,
    1525      e            t_seri,q_seri,u_seri,v_seri,soil_model,
    1526      e            ftsol,soilcap,soilflux,paprs,pplay,radsol,
    1527      e            fsnow,fqsol,
    1528      e            rlat, frugs,
     1555     e            t_seri,q_seri,u_seri,v_seri,ok_veget,
     1556     e            ftsol,paprs,pplay,radsol,
     1557     e            fsnow,fqsol,fevap,falbe,
     1558     e            rain_fall, snow_fall, solsw, sollw,
     1559     e            rlon, rlat, frugs,
     1560     e            debut, lafin,
    15291561     s            d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
    15301562     s            fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,rugmer,
     
    15321564     s            ycoefh,yu1,yv1)
    15331565c
    1534       DO i = 1, klon
    1535          sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
    1536          evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
     1566C§§§ PB
     1567C§§§ Incrementation des flux
     1568C§§
     1569      zxfluxt=0.
     1570      zxfluxq=0.
     1571      zxfluxu=0.
     1572      zxfluxv=0.
     1573      DO nsrf = 1, nbsrf
     1574        DO k = 1, klev
     1575          DO i = 1, klon
     1576            zxfluxt(i,k) = zxfluxt(i,k) +
     1577     $          fluxt(i,k,nsrf) * pctsrf( i, nsrf)
     1578            zxfluxq(i,k) = zxfluxq(i,k) +
     1579     $          fluxq(i,k,nsrf) * pctsrf( i, nsrf)
     1580            zxfluxu(i,k) = zxfluxu(i,k) +
     1581     $          fluxu(i,k,nsrf) * pctsrf( i, nsrf)
     1582            zxfluxv(i,k) = zxfluxv(i,k) +
     1583     $          fluxv(i,k,nsrf) * pctsrf( i, nsrf)
     1584          END DO
     1585        END DO
     1586      END DO
     1587      DO i = 1, klon
     1588         sens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
     1589c         evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
     1590         evap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
    15371591         fder(i) = dsens(i) + devap(i)
    15381592      ENDDO
     
    15511605      DO i = 1, klon
    15521606         zxtsol(i) = 0.0
     1607         IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
     1608     $       pctsrf(i, is_oce) + pctsrf(i, is_sic)  - 1.) .GT. EPSFRA)
     1609     $       THEN
     1610             WRITE(*,*) 'physiq : pb sous surface au point ', i,
     1611     $           pctsrf(i, 1 : nbsrf)
     1612         ENDIF
    15531613      ENDDO
    15541614      DO nsrf = 1, nbsrf
     
    15681628      ENDDO
    15691629
    1570 c
    1571 c Appeler le modele du sol
    1572 c
    1573       IF (soil_model) THEN
    1574          DO nsrf = 1, nbsrf
    1575          CALL soil(dtime, nsrf, fsnow(1,nsrf),
    1576      .             ftsol(1,nsrf), ftsoil(1,1,nsrf),
    1577      .             soilcap(1,nsrf), soilflux(1,nsrf))
    1578          ENDDO
    1579       ENDIF
    15801630c
    15811631c Calculer la derive du flux infrarouge
     
    16231673      ELSE IF (iflag_con.EQ.2) THEN
    16241674      CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
    1625      e            conv_t, conv_q, fluxq(1,1), omega,
     1675     e            conv_t, conv_q, zxfluxq(1,1), omega,
    16261676     s            d_t_con, d_q_con, rain_con, snow_con,
    16271677     s            pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
     
    18201870      CALL albsno(veget,agesno,alb_neig)
    18211871      DO i = 1, klon
    1822          zx_alb_oce = alb_eau(i)
     1872         falbe(i,is_oce) = alb_eau(i)
    18231873         IF (pctsrf(i,is_oce).GT.epsfra .AND. ftsol(i,is_oce).LT.271.35)
    1824      .   zx_alb_oce = 0.6 ! pour slab_ocean
     1874     .   falbe(i,is_oce) = 0.6 ! pour slab_ocean
    18251875         zfra = MAX(0.0,MIN(1.0,fsnow(i,is_lic)/(fsnow(i,is_lic)+10.0)))
    1826          zx_alb_lic = alb_neig(i)*zfra + 0.6*(1.0-zfra)
     1876         falbe(i,is_lic) = alb_neig(i)*zfra + 0.6*(1.0-zfra)
    18271877         zfra = MAX(0.0,MIN(1.0,fsnow(i,is_ter)/(fsnow(i,is_ter)+10.0)))
    1828          zx_alb_ter = alb_neig(i)*zfra + lmt_alb(i)*(1.0-zfra)
     1878         falbe(i,is_ter) = alb_neig(i)*zfra + lmt_alb(i)*(1.0-zfra)
    18291879         zfra = MAX(0.0,MIN(1.0,fsnow(i,is_sic)/(fsnow(i,is_sic)+10.0)))
    1830          zx_alb_sic = alb_neig(i)*zfra + 0.6*(1.0-zfra)
    1831          albsol(i) = zx_alb_oce * pctsrf(i,is_oce)
    1832      .             + zx_alb_lic * pctsrf(i,is_lic)
    1833      .             + zx_alb_ter * pctsrf(i,is_ter)
    1834      .             + zx_alb_sic * pctsrf(i,is_sic)
    1835       ENDDO
     1880         falbe(i,is_sic) = alb_neig(i)*zfra + 0.6*(1.0-zfra)
     1881         albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
     1882     .             + falbe(i,is_lic) * pctsrf(i,is_lic)
     1883     .             + falbe(i,is_ter) * pctsrf(i,is_ter)
     1884     .             + falbe(i,is_sic) * pctsrf(i,is_sic)
     1885      ENDDO
     1886c      DO nsrf = 1, nbsrf
     1887c        DO i = 1, klon
     1888c           albsol(i) = albsol(i) + falbe(i,nsrf)*pctsrf(i,nsrf)
     1889c        ENDDO
     1890c      ENDDO
    18361891      CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
    18371892     e            (dist, rmu0, fract, co2_ppm, solaire,
     
    18561911c Calculer l'hydrologie de la surface
    18571912c
    1858       CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, evap,
    1859      .            agesno, ftsol,fqsol,fsnow, ruis)
     1913c      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
     1914c     .            agesno, ftsol,fqsol,fsnow, ruis)
    18601915c
    18611916      DO i = 1, klon
     
    20072062c Accumuler les variables a stocker dans les fichiers histoire:
    20082063c
    2009       IF (ok_oasis) THEN ! couplage oasis
    2010       DO i = 1, klon
    2011         oas_sols(i) = oas_sols(i) + solsw(i)          / FLOAT(nexca)
    2012         oas_nsol(i) = oas_nsol(i) + (bils(i)-solsw(i))/ FLOAT(nexca)
    2013         oas_rain(i) = oas_rain(i) + rain_fall(i)      / FLOAT(nexca)
    2014         oas_snow(i) = oas_snow(i) + snow_fall(i)      / FLOAT(nexca)
    2015         oas_evap(i) = oas_evap(i) + evap(i)           / FLOAT(nexca)
    2016         oas_tsol(i) = oas_tsol(i) + zxtsol(i)         / FLOAT(nexca)
    2017         oas_fder(i) = oas_fder(i) + fder(i)           / FLOAT(nexca)
    2018         oas_albe(i) = oas_albe(i) + albsol(i)         / FLOAT(nexca)
    2019         oas_taux(i) = oas_taux(i) + fluxu(i,1)        / FLOAT(nexca)
    2020         oas_tauy(i) = oas_tauy(i) + fluxv(i,1)        / FLOAT(nexca)
    2021         oas_ruis(i) = oas_ruis(i) + ruis(i)       /FLOAT(nexca)/dtime
    2022       ENDDO
    2023       ENDIF
    20242064c
    20252065c
     
    20432083      CALL histwrite(nid_day,"tsol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
    20442084c
     2085C
     2086      zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, is_ter)
     2087      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d ,zx_tmp_2d)
     2088      CALL histwrite(nid_day,"tter",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
     2089C
     2090      zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, is_lic)
     2091      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
     2092      CALL histwrite(nid_day,"tlic",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
     2093C
     2094      zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, is_oce)
     2095      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
     2096      CALL histwrite(nid_day,"toce",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
     2097C
     2098      zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, is_sic)
     2099      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
     2100      CALL histwrite(nid_day,"tsic",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
     2101C
    20452102      DO i = 1, klon
    20462103         zx_tmp_fi2d(i) = paprs(i,1)
     
    20852142      CALL histwrite(nid_day,"ruis",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
    20862143c
    2087       DO i = 1, klon
    2088          zx_tmp_fi2d(i) = fluxu(i,1)
    2089       ENDDO
    2090       CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
    2091       CALL histwrite(nid_day,"frtu",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
    2092 c
    2093       DO i = 1, klon
    2094          zx_tmp_fi2d(i) = fluxv(i,1)
    2095       ENDDO
    2096       CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
    2097       CALL histwrite(nid_day,"frtv",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
    2098 c
    2099       DO i = 1, klon
    2100          zx_tmp_fi2d(i) = pctsrf(i,is_sic)
    2101       ENDDO
    2102       CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
    2103       CALL histwrite(nid_day,"sicf",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
     2144c      DO i = 1, klon
     2145c         zx_tmp_fi2d(i) = fluxu(i,1)
     2146c      ENDDO
     2147c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
     2148c      CALL histwrite(nid_day,"frtu",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
     2149c
     2150c      DO i = 1, klon
     2151c         zx_tmp_fi2d(i) = fluxv(i,1)
     2152c      ENDDO
     2153c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
     2154c      CALL histwrite(nid_day,"frtv",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
     2155c
     2156      DO nsrf = 1, nbsrf
     2157C§§§
     2158        zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)
     2159        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
     2160        CALL histwrite(nid_day,"pourc_"//clnsurf(nsrf),itap,
     2161     $      zx_tmp_2d,iim*jjmp1,ndex2d)
     2162C
     2163        zx_tmp_fi2d(1 : klon) = - fluxt( 1 : klon, 1, nsrf)
     2164        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
     2165        CALL histwrite(nid_day,"sens_"//clnsurf(nsrf),itap,
     2166     $      zx_tmp_2d,iim*jjmp1,ndex2d)
     2167C
     2168        zx_tmp_fi2d(1 : klon) = - fluxq( 1 : klon, 1, nsrf)
     2169        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
     2170        CALL histwrite(nid_day,"lat_"//clnsurf(nsrf),itap,
     2171     $      zx_tmp_2d,iim*jjmp1,ndex2d)
     2172C
     2173        zx_tmp_fi2d(1 : klon) = - fluxu( 1 : klon, 1, nsrf)
     2174        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
     2175        CALL histwrite(nid_day,"taux_"//clnsurf(nsrf),itap,
     2176     $      zx_tmp_2d,iim*jjmp1,ndex2d)
     2177C     
     2178        zx_tmp_fi2d(1 : klon) = - fluxv( 1 : klon, 1, nsrf)
     2179        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
     2180        CALL histwrite(nid_day,"tauy_"//clnsurf(nsrf),itap,
     2181     $      zx_tmp_2d,iim*jjmp1,ndex2d)
     2182C
     2183      END DO 
     2184C
     2185c$$$      DO i = 1, klon
     2186c$$$         zx_tmp_fi2d(i) = pctsrf(i,is_sic)
     2187c$$$      ENDDO
     2188c$$$      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
     2189c$$$      CALL histwrite(nid_day,"sicf",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
    21042190c
    21052191      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldl,zx_tmp_2d)
     
    22432329      CALL histwrite(nid_mth,"ruis",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
    22442330c
    2245       DO i = 1, klon
    2246          zx_tmp_fi2d(i) = fluxu(i,1)
    2247       ENDDO
    2248       CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
    2249       CALL histwrite(nid_mth,"frtu",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
    2250 c
    2251       DO i = 1, klon
    2252          zx_tmp_fi2d(i) = fluxv(i,1)
    2253       ENDDO
    2254       CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
    2255       CALL histwrite(nid_mth,"frtv",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
    2256 c
    2257       DO i = 1, klon
    2258          zx_tmp_fi2d(i) = pctsrf(i,is_sic)
    2259       ENDDO
    2260       CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
    2261       CALL histwrite(nid_mth,"sicf",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
     2331c      DO i = 1, klon
     2332c         zx_tmp_fi2d(i) = fluxu(i,1)
     2333c      ENDDO
     2334c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
     2335c      CALL histwrite(nid_mth,"frtu",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
     2336c
     2337c      DO i = 1, klon
     2338c         zx_tmp_fi2d(i) = fluxv(i,1)
     2339c      ENDDO
     2340c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
     2341c      CALL histwrite(nid_mth,"frtv",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
     2342c
     2343      DO nsrf = 1, nbsrf
     2344C§§§
     2345        zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)
     2346        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
     2347        CALL histwrite(nid_mth,"pourc_"//clnsurf(nsrf),itap,
     2348     $      zx_tmp_2d,iim*jjmp1,ndex2d)
     2349C
     2350        zx_tmp_fi2d(1 : klon) = - fluxt( 1 : klon, 1, nsrf)
     2351        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
     2352        CALL histwrite(nid_mth,"sens_"//clnsurf(nsrf),itap,
     2353     $      zx_tmp_2d,iim*jjmp1,ndex2d)
     2354C
     2355        zx_tmp_fi2d(1 : klon) = - fluxq( 1 : klon, 1, nsrf)
     2356        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
     2357        CALL histwrite(nid_mth,"lat_"//clnsurf(nsrf),itap,
     2358     $      zx_tmp_2d,iim*jjmp1,ndex2d)
     2359C
     2360        zx_tmp_fi2d(1 : klon) = - fluxu( 1 : klon, 1, nsrf)
     2361        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
     2362        CALL histwrite(nid_mth,"taux_"//clnsurf(nsrf),itap,
     2363     $      zx_tmp_2d,iim*jjmp1,ndex2d)
     2364C     
     2365        zx_tmp_fi2d(1 : klon) = - fluxv( 1 : klon, 1, nsrf)
     2366        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
     2367        CALL histwrite(nid_mth,"tauy_"//clnsurf(nsrf),itap,
     2368     $      zx_tmp_2d,iim*jjmp1,ndex2d)
     2369C
     2370      END DO 
     2371c$$$      DO i = 1, klon
     2372c$$$         zx_tmp_fi2d(i) = pctsrf(i,is_sic)
     2373c$$$      ENDDO
     2374c$$$      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
     2375c$$$      CALL histwrite(nid_mth,"sicf",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
    22622376c
    22632377      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsol,zx_tmp_2d)
     
    24762590      CALL histwrite(nid_ins,"psol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
    24772591c
     2592      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
     2593      CALL histwrite(nid_ins,"tsol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
     2594c
    24782595      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw,zx_tmp_2d)
    24792596      CALL histwrite(nid_ins,"topl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
     
    25082625      CALL gr_fi_ecrit(1, klon,iim,jjmp1, d_ts(1,is_sic),zx_tmp_2d)
    25092626      CALL histwrite(nid_ins,"dtsvdfi",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
     2627
     2628      DO nsrf = 1, nbsrf
     2629C§§§
     2630        zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)
     2631        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
     2632        CALL histwrite(nid_ins,"pourc_"//clnsurf(nsrf),itap,
     2633     $      zx_tmp_2d,iim*jjmp1,ndex2d)
     2634C
     2635        zx_tmp_fi2d(1 : klon) = - fluxt( 1 : klon, 1, nsrf)
     2636        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
     2637        CALL histwrite(nid_ins,"sens_"//clnsurf(nsrf),itap,
     2638     $      zx_tmp_2d,iim*jjmp1,ndex2d)
     2639C
     2640        zx_tmp_fi2d(1 : klon) = - fluxq( 1 : klon, 1, nsrf)
     2641        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
     2642        CALL histwrite(nid_ins,"lat_"//clnsurf(nsrf),itap,
     2643     $      zx_tmp_2d,iim*jjmp1,ndex2d)
     2644C
     2645        zx_tmp_fi2d(1 : klon) = ftsol( 1 : klon, nsrf)
     2646        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
     2647        CALL histwrite(nid_ins,"tsol_"//clnsurf(nsrf),itap,
     2648     $      zx_tmp_2d,iim*jjmp1,ndex2d)
     2649C
     2650        zx_tmp_fi2d(1 : klon) = - fluxu( 1 : klon, 1, nsrf)
     2651        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
     2652        CALL histwrite(nid_ins,"taux_"//clnsurf(nsrf),itap,
     2653     $      zx_tmp_2d,iim*jjmp1,ndex2d)
     2654C     
     2655        zx_tmp_fi2d(1 : klon) = - fluxv( 1 : klon, 1, nsrf)
     2656        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
     2657        CALL histwrite(nid_ins,"tauy_"//clnsurf(nsrf),itap,
     2658     $      zx_tmp_2d,iim*jjmp1,ndex2d)
     2659C
     2660      END DO 
    25102661
    25112662c
     
    25462697      ENDIF
    25472698c
    2548       IF (ok_oasis .AND. mod(itap,nexca).EQ.0) THEN
    2549 c
    2550 c Je ne traite pas le ruissellement, pour l'instant (qui m'aidera ?)
    2551          DO i = 1, klon
    2552             oas_ruisoce(i) = 0.0
    2553             oas_ruisriv(i) = 0.0
    2554          ENDDO
    2555 c
    2556          ig = 1
    2557          DO i = 1, iim
    2558             z_sols(i,1) = oas_sols(ig)
    2559             z_nsol(i,1) = oas_nsol(ig)
    2560             z_rain(i,1) = oas_rain(ig)
    2561             z_snow(i,1) = oas_snow(ig)
    2562             z_evap(i,1) = oas_evap(ig)
    2563             z_ruisoce(i,1) = oas_ruisoce(ig)
    2564             z_ruisriv(i,1) = oas_ruisriv(ig)
    2565             z_tsol(i,1) = oas_tsol(ig)
    2566             z_fder(i,1) = oas_fder(ig)
    2567             z_albe(i,1) = oas_albe(ig)
    2568             z_taux(i,1) = oas_taux(ig)
    2569             z_tauy(i,1) = oas_tauy(ig)
    2570          ENDDO
    2571          DO j = 2, jjm
    2572          DO i = 1, iim
    2573             ig = ig + 1
    2574             z_sols(i,j) = oas_sols(ig)
    2575             z_nsol(i,j) = oas_nsol(ig)
    2576             z_rain(i,j) = oas_rain(ig)
    2577             z_snow(i,j) = oas_snow(ig)
    2578             z_evap(i,j) = oas_evap(ig)
    2579             z_ruisoce(i,j) = oas_ruisoce(ig)
    2580             z_ruisriv(i,j) = oas_ruisriv(ig)
    2581             z_tsol(i,j) = oas_tsol(ig)
    2582             z_fder(i,j) = oas_fder(ig)
    2583             z_albe(i,j) = oas_albe(ig)
    2584             z_taux(i,j) = oas_taux(ig)
    2585             z_tauy(i,j) = oas_tauy(ig)
    2586          ENDDO
    2587          ENDDO
    2588          ig = ig + 1
    2589          DO i = 1, iim
    2590             z_sols(i,jjmp1)    = oas_sols(ig)
    2591             z_nsol(i,jjmp1)    = oas_nsol(ig)
    2592             z_rain(i,jjmp1)    = oas_rain(ig)
    2593             z_snow(i,jjmp1)    = oas_snow(ig)
    2594             z_evap(i,jjmp1)    = oas_evap(ig)
    2595             z_ruisoce(i,jjmp1) = oas_ruisoce(ig)
    2596             z_ruisriv(i,jjmp1) = oas_ruisriv(ig)
    2597             z_tsol(i,jjmp1)    = oas_tsol(ig)
    2598             z_fder(i,jjmp1)    = oas_fder(ig)
    2599             z_albe(i,jjmp1)    = oas_albe(ig)
    2600             z_taux(i,jjmp1)    = oas_taux(ig)
    2601             z_tauy(i,jjmp1)    = oas_tauy(ig)
    2602          ENDDO
    2603 c
    2604 c Passer les champs au coupleur:
    2605 c
    2606          CALL intocpl(itap,jjmp1*iim,
    2607      .                   z_sols, z_nsol,
    2608      .                   z_rain, z_snow, z_evap,
    2609      .                   z_ruisoce, z_ruisriv,
    2610      .                   z_tsol, z_fder, z_albe,
    2611      .                   z_taux, z_tauy)
    2612          DO i = 1, klon
    2613            oas_sols(i) = 0.0
    2614            oas_nsol(i) = 0.0
    2615            oas_rain(i) = 0.0
    2616            oas_snow(i) = 0.0
    2617            oas_evap(i) = 0.0
    2618            oas_ruis(i) = 0.0
    2619            oas_tsol(i) = 0.0
    2620            oas_fder(i) = 0.0
    2621            oas_albe(i) = 0.0
    2622            oas_taux(i) = 0.0
    2623            oas_tauy(i) = 0.0
    2624          ENDDO
    2625       ENDIF
    26262699c
    26272700c Ecrire la bande regionale (binaire grads)
     
    26392712         CALL ecriregs(84,bils)
    26402713         CALL ecriregs(84,pctsrf(1,is_sic))
    2641          CALL ecriregs(84,fluxu(1,1))
    2642          CALL ecriregs(84,fluxv(1,1))
     2714         CALL ecriregs(84,zxfluxu(1,1))
     2715         CALL ecriregs(84,zxfluxv(1,1))
    26432716         CALL ecriregs(84,ue)
    26442717         CALL ecriregs(84,ve)
     
    27052778ccc         IF (ok_oasis) CALL quitcpl
    27062779         CALL phyredem ("restartphy.nc",dtime,radpas,co2_ppm,solaire,
    2707      .        rlat,rlon,ftsol,ftsoil,deltat,fqsol,fsnow,
    2708      .        radsol,rugmer,agesno,
    2709      .        zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
    2710      .        t_ancien, q_ancien)
     2780     .      rlat, rlon, pctsrf, ftsol, ftsoil, deltat, fqsol, fsnow,
     2781     .      falbe, fevap, rain_fall, snow_fall,
     2782     .      solsw, sollw,
     2783     .      radsol,rugmer,agesno,
     2784     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
     2785     .      t_ancien, q_ancien)
    27112786      ENDIF
    27122787
Note: See TracChangeset for help on using the changeset viewer.