Ignore:
Timestamp:
Sep 4, 2015, 6:48:31 PM (9 years ago)
Author:
oboucher
Message:

Recoded diurnal cycle in physiq.F90
iflag_diurnal_cycle is the new flag
=0 for no diurnal cycle
=1 for old diurnal cycle
=2 for new diurnal cycle with updates
in between two calls to radiation according
to how rmu0 has changed over the previous
timestep of physiq

Location:
LMDZ5/trunk/libf/phylmd
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/orbite.F90

    r2346 r2359  
    102102END SUBROUTINE angle
    103103! ====================================================================
    104 SUBROUTINE zenang(longi, gmtime, pdtrad, lat, long, pmu0, frac)
     104SUBROUTINE zenang(longi, gmtime, pdtrad1, pdtrad2, lat, long, pmu0, frac)
    105105  USE dimphy
    106106  IMPLICIT NONE
    … …  
    114114  ! fournit des moyennes de pmu0 et non des valeurs
    115115  ! instantanees, du coup frac prend toutes les valeurs
    116   ! entre 0 et 1.
     116  ! entre 0 et 1. La routine integre entre gmtime+pdtrad1 et
     117  ! gmtime+pdtrad2 avec pdtrad1 et pdtrad2 exprimes en secondes.
    117118  ! Date   : premiere version le 13 decembre 1994
    118119  ! revu pour  GCM  le 30 septembre 1996
     120  ! revu le 3 septembre 2015 pour les bornes de l'integrale
    119121  ! ===============================================================
    120122  ! longi : la longitude vraie de la terre dans son plan
    121123  ! solaire a partir de l'equinoxe de printemps (degre)
    122124  ! gmtime : temps universel en fraction de jour
    123   ! pdtrad : pas de temps du rayonnement (secondes)
     125  ! pdtrad1 : borne inferieure du pas de temps du rayonnement (secondes)
     126  ! pdtrad2 : borne inferieure du pas de temps du rayonnement (secondes)
     127  ! pdtrad2-pdtrad1 correspond a pdtrad, le pas de temps du rayonnement (secondes)
    124128  ! lat------INPUT : latitude en degres
    125129  ! long-----INPUT : longitude en degres
    126   ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
    127   ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
    128   ! ================================================================
    129   include "YOMCST.h"
    130   ! ================================================================
    131   REAL, INTENT (IN) :: longi, gmtime, pdtrad
     130  ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime+pdtrad1 et gmtime+pdtrad2
     131  ! frac-----OUTPUT: ensoleillement moyen entre gmtime+pdtrad1 et gmtime+pdtrad2
     132  ! ================================================================
     133  include "YOMCST.h"
     134  ! ================================================================
     135  REAL, INTENT (IN) :: longi, gmtime, pdtrad1, pdtrad2
    132136  REAL lat(klon), long(klon), pmu0(klon), frac(klon)
    133137  ! ================================================================
    … …  
    143147  REAL lat_sun ! declinaison en radian
    144148  REAL lon_sun ! longitude solaire en radian
    145   REAL latr ! latitude du pt de grille en radian
     149  REAL latr    ! latitude du pt de grille en radian
    146150  ! ================================================================
    147151
    … …  
    153157  lat_sun = asin(sin(lon_sun)*sin(incl))
    154158
    155   gmtime1 = gmtime*86400.
    156   gmtime2 = gmtime*86400. + pdtrad
     159  gmtime1 = gmtime*86400. + pdtrad1
     160  gmtime2 = gmtime*86400. + pdtrad2
    157161
    158162  DO i = 1, klon
    … …  
    160164    latr = lat(i)*pi_local/180.
    161165
    162     ! --pose probleme quand lat=+/-90 degres
    163 
    164     ! omega = -TAN(latr)*TAN(lat_sun)
    165     ! omega = ACOS(omega)
    166     ! IF (latr.GE.(pi_local/2.+lat_sun)
    167     ! .    .OR. latr.LE.(-pi_local/2.+lat_sun)) THEN
    168     ! omega = 0.0       ! nuit polaire
    169     ! ENDIF
    170     ! IF (latr.GE.(pi_local/2.-lat_sun)
    171     ! .          .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN
    172     ! omega = pi_local  ! journee polaire
    173     ! ENDIF
    174 
    175     ! --remplace par cela (le cas par defaut est different)
    176 
    177166    omega = 0.0 !--nuit polaire
     167
    178168    IF (latr>=(pi_local/2.-lat_sun) .OR. latr<=(-pi_local/2.-lat_sun)) THEN
    179169      omega = pi_local ! journee polaire
    180170    END IF
     171
    181172    IF (latr<(pi_local/2.+lat_sun) .AND. latr>(-pi_local/2.+lat_sun) .AND. &
    182173        latr<(pi_local/2.-lat_sun) .AND. latr>(-pi_local/2.-lat_sun)) THEN

  • LMDZ5/trunk/libf/phylmd/phys_state_var_mod.F90

    r2345 r2359  
    9797!$OMP THREADPRIVATE(clwcon0th,rnebcon0th)
    9898! radiation outputs
     99      REAL,ALLOCATABLE,SAVE :: zswdn0(:,:), zswdn(:,:)
     100!$OMP THREADPRIVATE(zswdn0,zswdn)
     101      REAL,ALLOCATABLE,SAVE :: zswup0(:,:), zswup(:,:)
     102!$OMP THREADPRIVATE(zswup0,zswup)
    99103      REAL,ALLOCATABLE,SAVE :: swdn0(:,:), swdn(:,:)
    100104!$OMP THREADPRIVATE(swdn0,swdn)
    … …  
    340344      REAL,ALLOCATABLE,SAVE :: lwup0p(:,:), lwupp(:,:)
    341345!$OMP THREADPRIVATE(lwdn0p, lwdnp, lwup0p, lwupp)
     346      REAL,ALLOCATABLE,SAVE :: zswdn0p(:,:), zswdnp(:,:)
     347      REAL,ALLOCATABLE,SAVE :: zswup0p(:,:), zswupp(:,:)
     348!$OMP THREADPRIVATE(zswdn0p, zswdnp, zswup0p, zswupp)
    342349      REAL,ALLOCATABLE,SAVE :: swdn0p(:,:), swdnp(:,:)
    343350      REAL,ALLOCATABLE,SAVE :: swup0p(:,:), swupp(:,:)
    … …  
    458465      ALLOCATE(clwcon0th(klon,klev),rnebcon0th(klon,klev))
    459466! radiation outputs
     467      ALLOCATE(zswdn0(klon,klevp1), zswdn(klon,klevp1))
     468      ALLOCATE(zswup0(klon,klevp1), zswup(klon,klevp1))
    460469      ALLOCATE(swdn0(klon,klevp1), swdn(klon,klevp1))
    461470      ALLOCATE(swup0(klon,klevp1), swup(klon,klevp1))
    … …  
    557566      ALLOCATE(lwdn0p(klon,klevp1), lwdnp(klon,klevp1))
    558567      ALLOCATE(lwup0p(klon,klevp1), lwupp(klon,klevp1))
     568      ALLOCATE(zswdn0p(klon,klevp1), zswdnp(klon,klevp1))
     569      ALLOCATE(zswup0p(klon,klevp1), zswupp(klon,klevp1))
    559570      ALLOCATE(swdn0p(klon,klevp1), swdnp(klon,klevp1))
    560571      ALLOCATE(swup0p(klon,klevp1), swupp(klon,klevp1))
    … …  
    611622      deallocate(clwcon0th, rnebcon0th)
    612623! radiation outputs
     624      deallocate(zswdn0, zswdn)
     625      deallocate(zswup0, zswup)
    613626      deallocate(swdn0, swdn)
    614627      deallocate(swup0, swup)
    … …  
    687700      deallocate(lwdn0p, lwdnp)
    688701      deallocate(lwup0p, lwupp)
     702      deallocate(zswdn0p, zswdnp)
     703      deallocate(zswup0p, zswupp)
    689704      deallocate(swdn0p, swdnp)
    690705      deallocate(swup0p, swupp)

  • LMDZ5/trunk/libf/phylmd/physiq.F90

    r2357 r2359  
    143143  !======================================================================
    144144  ! Clef controlant l'activation du cycle diurne:
    145   !cc      LOGICAL cycle_diurne
    146   !cc      PARAMETER (cycle_diurne=.FALSE.)
     145  ! en attente du codage des cles par Fred
     146        INTEGER iflag_cycle_diurne
     147        PARAMETER (iflag_cycle_diurne=1)
    147148  !======================================================================
    148149  ! Modele thermique du sol, a activer pour le cycle diurne:
    … …  
    589590  !
    590591  REAL dist, rmu0(klon), fract(klon)
    591   REAL zdtime, zlongi
     592  REAL zrmu0(klon), zfract(klon), swradcorr(klon)
     593  REAL zdtime, zdtime1, zdtime2, zlongi
    592594  !
    593595  REAL qcheck
    … …  
    10701072     print*,'iflag_coupl,iflag_clos,iflag_wake', &
    10711073          iflag_coupl,iflag_clos,iflag_wake
    1072      print*,'CYCLE_DIURNE', cycle_diurne
     1074     print*,'iflag_CYCLE_DIURNE', iflag_cycle_diurne
    10731075     !
    10741076     IF (iflag_con.EQ.2.AND.iflag_cld_th.GT.-1) THEN
    … …  
    11221124     ENDIF
    11231125
    1124      !IM cf. AM 081204 BEG
    1125      PRINT*,'cycle_diurne3 =',cycle_diurne
    1126      !IM cf. AM 081204 END
    1127      !
    11281126     CALL printflag( tabcntr0,radpas,ok_journe, &
    11291127          ok_instan, ok_region )
    … …  
    11491147     ENDIF
    11501148     !
    1151      IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
     1149     IF (dtime*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN
    11521150        WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
    11531151        WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
    … …  
    18211819  ! non nul aux poles.
    18221820  IF (abs(solarlong0-1000.)<1.e-4) then
    1823      call zenang_an(cycle_diurne,jH_cur,rlat,rlon,rmu0,fract)
     1821     call zenang_an(iflag_cycle_diurne.GE.1,jH_cur,rlat,rlon,rmu0,fract)
    18241822  ELSE
    1825      !  Avec ou sans cycle diurne
    1826      IF (cycle_diurne) THEN
     1823  ! recode par Olivier Boucher en sept 2015
     1824     SELECT CASE (iflag_cycle_diurne)
     1825     CASE(0) 
     1826     !  Sans cycle diurne
     1827        CALL angle(zlongi, rlat, fract, rmu0)
     1828     CASE(1) 
     1829     !  Avec cycle diurne sans application des poids
     1830     !  bit comparable a l ancienne formulation cycle_diurne=true
     1831     !  on integre entre gmtime et gmtime+radpas
    18271832        zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
    1828         CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
    1829      ELSE
    1830         CALL angle(zlongi, rlat, fract, rmu0)
    1831      ENDIF
     1833        CALL zenang(zlongi,jH_cur,0.0,zdtime,rlat,rlon,rmu0,fract)
     1834     CASE(2) 
     1835     !  Avec cycle diurne sans application des poids
     1836     !  On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
     1837     !  Comme cette routine est appele a tous les pas de temps de la physique
     1838     !  meme si le rayonnement n'est pas appele je remonte en arriere les
     1839     !  radpas-1 pas de temps suivant. Petite ruse pour prendre en compte le
     1840     !  premier pas de temps la physique ou itaprad=0
     1841        zdtime1=dtime*REAL(-MOD(itaprad,4)-1)     
     1842        zdtime2=dtime*REAL(radpas-MOD(itaprad,4)-1)
     1843        CALL zenang(zlongi,jH_cur,zdtime1,zdtime2,rlat,rlon,rmu0,fract)
     1844     END SELECT
    18321845  ENDIF
    1833 
    1834   ! AI Janv 2014
    1835   do i = 1, klon
    1836      if (fract(i).le.0.) then
    1837         JrNt(i)=0.
    1838      else
    1839         JrNt(i)=1.
    1840      endif
    1841   enddo
    18421846
    18431847  if (mydebug) then
    … …  
    35623566             topsw0,toplw0,solsw0,sollw0, &
    35633567             lwdn0, lwdn, lwup0, lwup,  &
    3564              swdn0, swdn, swup0, swup, &
     3568             zswdn0, zswdn, zswup0, zswup, &
    35653569             ok_ade, ok_aie, &
    35663570             tau_aero, piz_aero, cg_aero, &
    … …  
    36083612             topsw0,toplw0,solsw0,sollw0, &
    36093613             lwdn0, lwdn, lwup0, lwup,  &
    3610              swdn0, swdn, swup0, swup, &
     3614             zswdn0, zswdn, zswup0, zswup, &
    36113615             topswad_aero, solswad_aero, &
    36123616             topswai_aero, solswai_aero, &
    … …  
    36643668                   topsw0p,toplw0p,solsw0p,sollw0p, &
    36653669                   lwdn0p, lwdnp, lwup0p, lwupp,  &
    3666                    swdn0p, swdnp, swup0p, swupp, &
     3670                   zswdn0p, zswdnp, zswup0p, zswupp, &
    36673671                   topswad_aerop, solswad_aerop, &
    36683672                   topswai_aerop, solswai_aerop, &
    … …  
    37063710     swup=0.    ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
    37073711     swup0=0.
    3708      swdn=0.
    3709      swdn0=0.
     3712     zswdn=0.
     3713     zswdn0=0.
    37103714     lwup=0.
    37113715     lwup0=0.
    … …  
    37153719
    37163720  !
     3721  ! Calculer les poids a appliquer sur le SW
     3722  ! sortie JrNt = jour-nuit
     3723  ! recode par Olivier Boucher en sept 2015
     3724  !
     3725
     3726  SELECT CASE (iflag_cycle_diurne)
     3727  !
     3728  CASE(0)
     3729  !  Sans cycle diurne
     3730     swradcorr=1.0
     3731     JrNt = 1.0
     3732  CASE(1)
     3733  !  Avec cycle diurne sans les poids
     3734     swradcorr=1.0
     3735     JrNt=0.0
     3736     WHERE (fract.GT.0.0) JrNt=1.0
     3737  CASE(2)
     3738  !  Avec cycle diurne et les poids
     3739     zdtime1=-dtime !--on corrige le rayonnement pour representer le
     3740     zdtime2=0.0    !--pas de temps de la physique qui se termine
     3741     CALL zenang(zlongi,jH_cur,zdtime1,zdtime2,rlat,rlon,zrmu0,zfract)
     3742     swradcorr=0.0
     3743     WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) swradcorr=zfract/fract*zrmu0/rmu0
     3744     JrNt=0.0
     3745     WHERE (zfract.GT.0.0) JrNt=1.0
     3746  END SELECT
     3747
     3748  !
     3749  ! Corriger les flux SW pour le cycle diurne ameliore
     3750  ! recode par Olivier Boucher en sept 2015
     3751  !
     3752
     3753  DO k=1, klev+1
     3754    swdn0(:,k)=swradcorr(:)*zswdn0(:,k)
     3755    swdn(:,k) =swradcorr(:)*zswdn(:,k)
     3756    swup0(:,k)=swradcorr(:)*zswup0(:,k)
     3757    swup(:,k) =swradcorr(:)*zswup(:,k)
     3758  ENDDO
     3759  if (ok_4xCO2atm) then
     3760  DO k=1, klev+1
     3761    swdn0p(:,k)=swradcorr(:)*zswdn0p(:,k)
     3762    swdnp(:,k) =swradcorr(:)*zswdnp(:,k)
     3763    swup0p(:,k)=swradcorr(:)*zswup0p(:,k)
     3764    swupp(:,k) =swradcorr(:)*zswupp(:,k)
     3765  ENDDO
     3766  endif
     3767
     3768  !
    37173769  ! Ajouter la tendance des rayonnements (tous les pas)
    37183770  !
    3719   d_t_swr(:,:)=heat(:,:)*dtime/RDAY
    3720   d_t_lwr(:,:)=-cool(:,:)*dtime/RDAY
    3721   d_t_sw0(:,:)=heat0(:,:)*dtime/RDAY
    3722   d_t_lw0(:,:)=-cool0(:,:)*dtime/RDAY
     3771 
     3772  DO k=1, klev
     3773    d_t_swr(:,k)=swradcorr(:)*heat(:,k)*dtime/RDAY
     3774    d_t_sw0(:,k)=swradcorr(:)*heat0(:,k)*dtime/RDAY
     3775    d_t_lwr(:,k)=-cool(:,k)*dtime/RDAY
     3776    d_t_lw0(:,k)=-cool0(:,k)*dtime/RDAY
     3777  ENDDO
     3778
    37233779  CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy)
    37243780  CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy)
Note: See TracChangeset for help on using the changeset viewer.