Ignore:
Timestamp:
Jul 7, 2009, 4:01:00 PM (15 years ago)
Author:
Laurent Fairhead
Message:

Modifications nécessaires a l'inclusion d'un calendrier réaliste.
La date courante est calculée dans leapfrog.F et exprimée en Jour Julien
(modifié). On en a profité pour faire un peu de ménage dans la gestion des dates
du modèle.
Dans la physique, on utilise les routines de passages entre calendrier Julien et
Gregorien incluses dans IOIPSL pour calculer le nombre de jours écoulés depuis le
1er janvier (pour les conditions aux limites) ou l'equinoxe (pour le calcul de
la longitude solaire). Le calcul de l'orbite reprend celui du gcm planétaire
(codé par FH)
On décide du calendrier à utiliser à l'aide du paramètre calend du run.def. Par
défaut celui-ci est à earth_360d
LF

Location:
LMDZ4/branches/LMDZ4-dev/libf/phylmd
Files:
2 added
2 edited

Legend:

Unmodified
Added
Removed
  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/hgardfou.F

    r1145 r1201  
    11!
     2! $Id$
    23      SUBROUTINE hgardfou (t,tsol,text)
    34      use dimphy
     
    1213      REAL t(klon,klev), tsol(klon,nbsrf)
    1314      CHARACTER*(*) text
     15      character (len=20) :: modname = 'hgardfou'
     16      character (len=80) :: abort_message
    1417C
    1518      INTEGER i, k, nsrf
     
    124127c
    125128      IF (.NOT. ok) THEN
    126          PRINT*, 'hgardfou s arrete ', text
    127          CALL abort
     129         abort_message= 'hgardfou s arrete '//text
     130         CALL abort_gcm (modname,abort_message,1)
    128131      ENDIF
    129132
  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F

    r1196 r1201  
    44
    55      SUBROUTINE physiq (nlon,nlev,
    6      .            debut,lafin,rjourvrai,gmtime,pdtphys,
     6     .            debut,lafin,jD_cur, jH_cur,pdtphys,
    77     .            paprs,pplay,pphi,pphis,presnivs,clesphy0,
    88     .            u,v,t,qx,
     
    6060c debut---input-L-variable logique indiquant le premier passage
    6161c lafin---input-L-variable logique indiquant le dernier passage
    62 c rjour---input-R-numero du jour de l'experience
    63 c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
     62c jD_cur--input-R-jour courant a l'appel de la physique (jour julien)
     63c jH_cur--input-R-heure courante a l'appel de la physique (jour julien)
    6464c pdtphys-input-R-pas d'integration pour la physique (seconde)
    6565c paprs---input-R-pression pour chaque inter-couche (en Pa)
     
    186186      INTEGER nlon
    187187      INTEGER nlev
    188       REAL rjourvrai
    189       REAL gmtime
     188      REAL :: jD_cur, jH_cur
     189
    190190      REAL pdtphys
    191191      LOGICAL debut, lafin
     
    703703c Conditions aux limites
    704704c
    705       INTEGER julien
     705!
     706! Gestion calendrier
     707!
     708      REAL :: jD_1jan, jH_1jan
     709      INTEGER :: year_cur, mth_cur, day_cur, days_elapsed
     710      REAL :: hour, day_since_equinox
     711! Date de l'equinoxe de printemps
     712      INTEGER, parameter :: mth_eq=3, day_eq=21
     713      REAL :: jD_eq
     714
     715      LOGICAL, parameter :: new_orbit = .true.
     716
    706717c
    707718      INTEGER lmt_pas
     
    10871098c$OMP THREADPRIVATE(first)
    10881099
     1100      integer iunit
     1101
    10891102      logical, save::  read_climoz ! read ozone climatology
    10901103      integer, save:: ncid_climoz ! NetCDF file containing ozone climatology
     
    11141127         write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
    11151128         write(lunout,*)
    1116      s 'nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys'
     1129     s 'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
    11171130         write(lunout,*)
    1118      s  nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys
     1131     s  nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
    11191132
    11201133         write(lunout,*) 'papers, play, phi, u, v, t, omega'
     
    11841197
    11851198c======================================================================
    1186       xjour = rjourvrai
     1199! Gestion calendrier
     1200!
     1201      call ju2ymds(jD_cur+jH_cur, year_cur, mth_cur, day_cur, hour)
     1202      call ymds2ju(year_cur, 1, 1, 0., jD_1jan)
     1203      jH_1jan = jD_1jan - int (jD_1jan)
     1204      jD_1jan = int (jD_1jan)
     1205      xjour = jD_cur - jD_1jan
     1206      days_elapsed = jD_cur - jD_1jan
     1207
    11871208c
    11881209c Si c'est le debut, il faut initialiser plusieurs choses
     
    14721493
    14731494cXXXPB Positionner date0 pour initialisation de ORCHIDEE
    1474       date0 = zjulian
    1475 C      date0 = day_ini
     1495      date0 = jD_ref
    14761496      WRITE(*,*) 'physiq date0 : ',date0
    14771497c
     
    14911511         CALL VTe(VTphysiq)
    14921512         CALL VTb(VTinca)
    1493          iii = MOD(NINT(xjour),360)
    1494          calday = FLOAT(iii) + gmtime
    1495          WRITE(lunout,*) 'initial time ', xjour, calday
     1513!         iii = MOD(NINT(xjour),360)
     1514!         calday = FLOAT(iii) + jH_cur
     1515         calday = FLOAT(days_elapsed) + jH_cur
     1516         WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
    14961517
    14971518         CALL chemini(
     
    15341555!
    15351556      itap   = itap + 1
    1536       julien = MOD(NINT(xjour),360)
    1537       if (julien .eq. 0) julien = 360
    1538 
    15391557!
    15401558! Update fraction of the sub-surfaces (pctsrf) and
     
    15421560! on the surface fraction.
    15431561!
    1544       CALL change_srf_frac(itap, dtime, julien,
     1562      CALL change_srf_frac(itap, dtime, days_elapsed+1,
    15451563     *     pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
    1546 
    15471564
    15481565! Tendances bidons pour les processus qui n'affectent pas certaines
     
    16971714         if (read_climoz) then
    16981715C           Ozone climatology from a NetCDF file
    1699             call regr_pr_av(ncid_climoz, "tro3", julien, press_climoz,
     1716            call regr_pr_av(ncid_climoz, "tro3", days_elapsed+1,
     1717     &           press_climoz,
    17001718     $           paprs, wo)
    17011719!           Convert from mole fraction of ozone to column density of ozone in a
     
    17071725C           "zmasse" changes a little.)
    17081726         else
    1709             CALL ozonecm(real(julien), rlat, paprs, wo)
     1727            CALL ozonecm(real(days_elapsed+1), rlat, paprs, wo)
     1728
    17101729         end if
    17111730      ENDIF
     
    17491768! doit donc etre placé avant radlwsw et pbl_surface
    17501769
     1770! calcul selon la routine utilisee pour les planetes
     1771      if (new_orbit) then
     1772        call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
     1773        day_since_equinox = (jD_cur + jH_cur) - jD_eq
     1774!        day_since_equinox = (jD_cur) - jD_eq
     1775        call solarlong(day_since_equinox, zlongi, dist)
     1776      else     
     1777! calcul selon la routine utilisee pour l'AR4
    17511778!   choix entre calcul de la longitude solaire vraie ou valeur fixee a
    17521779!   solarlong0
    1753 
    1754       if (solarlong0<-999.) then
    1755          CALL orbite(FLOAT(julien),zlongi,dist)
    1756       else
    1757          zlongi=solarlong0  ! longitude solaire vraie
    1758          dist=1.            ! distance au soleil / moyenne
     1780        if (solarlong0<-999.) then
     1781           CALL orbite(FLOAT(days_elapsed+1),zlongi,dist)
     1782        else
     1783           zlongi=solarlong0  ! longitude solaire vraie
     1784           dist=1.            ! distance au soleil / moyenne
     1785        endif
    17591786      endif
    1760 
    1761       if(prt_level.ge.1) print*,'Longitude solaire ',zlongi,solarlong0
     1787      if(prt_level.ge.1)                                                &
     1788     &    write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
    17621789
    17631790!  Avec ou sans cycle diurne
    17641791      IF (cycle_diurne) THEN
    17651792        zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
    1766         CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
     1793        CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
    17671794      ELSE
    17681795        CALL angle(zlongi, rlat, fract, rmu0)
     
    17971824
    17981825      CALL pbl_surface(
    1799      e     dtime,     date0,     itap,    julien,
     1826     e     dtime,     date0,     itap,    days_elapsed+1,
    18001827     e     debut,     lafin,
    18011828     e     rlon,      rlat,      rugoro,  rmu0,     
     
    20942121
    20952122            if (itop_con(i).gt.klev-3) then
    2096                print*,'La convection monte trop haut '
    2097                print*,'itop_con(,',i,',)=',itop_con(i)
     2123              if(prt_level >= 9) then
     2124                write(lunout,*)'La convection monte trop haut '
     2125                write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
     2126              endif
    20982127            endif
    20992128          ENDDO
     
    26352664         IF (.NOT. aerosol_couple)
    26362665     &        CALL readaerosol_optic(
    2637      &        debut, new_aod, flag_aerosol, rjourvrai, pdtphys,
     2666     &        debut, new_aod, flag_aerosol, jD_cur-jD_ref, pdtphys,
    26382667     &        pplay, paprs, t_seri, rhcl,
    26392668     &        mass_solu_aero, mass_solu_aero_pi,
     
    27642793         CALL VTe(VTphysiq)
    27652794         CALL VTb(VTinca)
    2766          calday = FLOAT(julien) + gmtime
     2795         calday = FLOAT(days_elapsed + 1) + jH_cur
    27672796
    27682797         IF (config_inca == 'aero') THEN
     
    27752804
    27762805         CALL chemhook_begin (calday,
    2777      $                          julien,
    2778      $                          gmtime,
     2806     $                          days_elapsed,
     2807     $                          jH_cur,
    27792808     $                          pctsrf(1,1),
    27802809     $                          rlat,
     
    31063135      IF (is_sequential) THEN
    31073136     
    3108         CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
     3137        CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
    31093138     C                 ra,rg,romega,
    31103139     C                 rlat,rlon,pphis,
     
    31333162
    31343163      call phytrac (
    3135      I     itap,     julien,    gmtime,   debut,
     3164     I     itap,     days_elapsed+1,    jH_cur,   debut,
    31363165     I     lafin,    dtime,     u, v,     t,
    31373166     I     paprs,    pplay,     pmfu,     pmfd,
     
    33513380      write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
    33523381      write(lunout,*)
    3353      s 'nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos'
     3382     s 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
    33543383      write(lunout,*)
    3355      s  nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys,
     3384     s  nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys,
    33563385     s  pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
    33573386     s  pctsrf(igout,is_sic)
     
    34473476      ENDIF
    34483477     
     3478!      first=.false.
    34493479
    34503480      RETURN
Note: See TracChangeset for help on using the changeset viewer.