Changeset 1201 for LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F

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

File:

• LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F (modified) (23 diffs)

LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F

@@ -4 +4 @@
 
       SUBROUTINE physiq (nlon,nlev,
-     .            debut,lafin,rjourvrai,gmtime,pdtphys,
+     .            debut,lafin,jD_cur, jH_cur,pdtphys,
      .            paprs,pplay,pphi,pphis,presnivs,clesphy0,
      .            u,v,t,qx,
@@ -60 +60 @@
 c debut---input-L-variable logique indiquant le premier passage
 c lafin---input-L-variable logique indiquant le dernier passage
-c rjour---input-R-numero du jour de l'experience
-c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
+c jD_cur--input-R-jour courant a l'appel de la physique (jour julien)
+c jH_cur--input-R-heure courante a l'appel de la physique (jour julien)
 c pdtphys-input-R-pas d'integration pour la physique (seconde)
 c paprs---input-R-pression pour chaque inter-couche (en Pa)
@@ -186 +186 @@
       INTEGER nlon
       INTEGER nlev
-      REAL rjourvrai
-      REAL gmtime
+      REAL :: jD_cur, jH_cur
+
       REAL pdtphys
       LOGICAL debut, lafin
@@ -703 +703 @@
 c Conditions aux limites
 c
-      INTEGER julien
+!
+! Gestion calendrier
+!
+      REAL :: jD_1jan, jH_1jan
+      INTEGER :: year_cur, mth_cur, day_cur, days_elapsed
+      REAL :: hour, day_since_equinox
+! Date de l'equinoxe de printemps
+      INTEGER, parameter :: mth_eq=3, day_eq=21
+      REAL :: jD_eq
+
+      LOGICAL, parameter :: new_orbit = .true.
+
 c
       INTEGER lmt_pas
@@ -1087 +1098 @@
 c$OMP THREADPRIVATE(first)
 
+      integer iunit
+
       logical, save::  read_climoz ! read ozone climatology
       integer, save:: ncid_climoz ! NetCDF file containing ozone climatology
@@ -1114 +1127 @@
          write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
          write(lunout,*)
-     s 'nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys'
+     s 'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
          write(lunout,*)
-     s  nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys 
+     s  nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys 
 
          write(lunout,*) 'papers, play, phi, u, v, t, omega'
@@ -1184 +1197 @@
 
 c======================================================================
-      xjour = rjourvrai
+! Gestion calendrier
+!
+      call ju2ymds(jD_cur+jH_cur, year_cur, mth_cur, day_cur, hour)
+      call ymds2ju(year_cur, 1, 1, 0., jD_1jan)
+      jH_1jan = jD_1jan - int (jD_1jan)
+      jD_1jan = int (jD_1jan) 
+      xjour = jD_cur - jD_1jan
+      days_elapsed = jD_cur - jD_1jan
+
 c
 c Si c'est le debut, il faut initialiser plusieurs choses
@@ -1472 +1493 @@
 
 cXXXPB Positionner date0 pour initialisation de ORCHIDEE
-      date0 = zjulian
-C      date0 = day_ini
+      date0 = jD_ref 
       WRITE(*,*) 'physiq date0 : ',date0
 c
@@ -1491 +1511 @@
          CALL VTe(VTphysiq)
          CALL VTb(VTinca)
-         iii = MOD(NINT(xjour),360)
-         calday = FLOAT(iii) + gmtime
-         WRITE(lunout,*) 'initial time ', xjour, calday
+!         iii = MOD(NINT(xjour),360)
+!         calday = FLOAT(iii) + jH_cur
+         calday = FLOAT(days_elapsed) + jH_cur
+         WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
 
          CALL chemini( 
@@ -1534 +1555 @@
 !
       itap   = itap + 1
-      julien = MOD(NINT(xjour),360)
-      if (julien .eq. 0) julien = 360
-
 !
 ! Update fraction of the sub-surfaces (pctsrf) and 
@@ -1542 +1560 @@
 ! on the surface fraction.
 !
-      CALL change_srf_frac(itap, dtime, julien, 
+      CALL change_srf_frac(itap, dtime, days_elapsed+1, 
      *     pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
-
 
 ! Tendances bidons pour les processus qui n'affectent pas certaines
@@ -1697 +1714 @@
          if (read_climoz) then
 C           Ozone climatology from a NetCDF file
-            call regr_pr_av(ncid_climoz, "tro3", julien, press_climoz,
+            call regr_pr_av(ncid_climoz, "tro3", days_elapsed+1, 
+     &           press_climoz,
      $           paprs, wo)
 !           Convert from mole fraction of ozone to column density of ozone in a
@@ -1707 +1725 @@
 C           "zmasse" changes a little.)
          else
-            CALL ozonecm(real(julien), rlat, paprs, wo)
+            CALL ozonecm(real(days_elapsed+1), rlat, paprs, wo)
+
          end if
       ENDIF
@@ -1749 +1768 @@
 ! doit donc etre placé avant radlwsw et pbl_surface
 
+! calcul selon la routine utilisee pour les planetes
+      if (new_orbit) then
+        call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
+        day_since_equinox = (jD_cur + jH_cur) - jD_eq 
+!        day_since_equinox = (jD_cur) - jD_eq 
+        call solarlong(day_since_equinox, zlongi, dist)
+      else      
+! calcul selon la routine utilisee pour l'AR4
 !   choix entre calcul de la longitude solaire vraie ou valeur fixee a
 !   solarlong0
-
-      if (solarlong0<-999.) then
-         CALL orbite(FLOAT(julien),zlongi,dist)
-      else
-         zlongi=solarlong0  ! longitude solaire vraie
-         dist=1.            ! distance au soleil / moyenne
+        if (solarlong0<-999.) then
+           CALL orbite(FLOAT(days_elapsed+1),zlongi,dist)
+        else
+           zlongi=solarlong0  ! longitude solaire vraie
+           dist=1.            ! distance au soleil / moyenne
+        endif
       endif
-
-      if(prt_level.ge.1) print*,'Longitude solaire ',zlongi,solarlong0
+      if(prt_level.ge.1)                                                &
+     &    write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
 
 !  Avec ou sans cycle diurne
       IF (cycle_diurne) THEN
         zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
-        CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
+        CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
       ELSE
         CALL angle(zlongi, rlat, fract, rmu0)
@@ -1797 +1824 @@
 
       CALL pbl_surface( 
-     e     dtime,     date0,     itap,    julien,
+     e     dtime,     date0,     itap,    days_elapsed+1,
      e     debut,     lafin,
      e     rlon,      rlat,      rugoro,  rmu0,     
@@ -2094 +2121 @@
 
             if (itop_con(i).gt.klev-3) then
-               print*,'La convection monte trop haut '
-               print*,'itop_con(,',i,',)=',itop_con(i)
+              if(prt_level >= 9) then
+                write(lunout,*)'La convection monte trop haut '
+                write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
+              endif
             endif
           ENDDO
@@ -2635 +2664 @@
          IF (.NOT. aerosol_couple)
      &        CALL readaerosol_optic(
-     &        debut, new_aod, flag_aerosol, rjourvrai, pdtphys,
+     &        debut, new_aod, flag_aerosol, jD_cur-jD_ref, pdtphys,
      &        pplay, paprs, t_seri, rhcl, 
      &        mass_solu_aero, mass_solu_aero_pi, 
@@ -2764 +2793 @@
          CALL VTe(VTphysiq)
          CALL VTb(VTinca)
-         calday = FLOAT(julien) + gmtime
+         calday = FLOAT(days_elapsed + 1) + jH_cur
 
          IF (config_inca == 'aero') THEN
@@ -2775 +2804 @@
 
          CALL chemhook_begin (calday,
-     $                          julien,
-     $                          gmtime,
+     $                          days_elapsed,
+     $                          jH_cur,
      $                          pctsrf(1,1),
      $                          rlat,
@@ -3106 +3135 @@
       IF (is_sequential) THEN
       
-        CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
+        CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
      C                 ra,rg,romega,
      C                 rlat,rlon,pphis,
@@ -3133 +3162 @@
 
       call phytrac (
-     I     itap,     julien,    gmtime,   debut,
+     I     itap,     days_elapsed+1,    jH_cur,   debut,
      I     lafin,    dtime,     u, v,     t,
      I     paprs,    pplay,     pmfu,     pmfd, 
@@ -3351 +3380 @@
       write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
       write(lunout,*)
-     s 'nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos'
+     s 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
       write(lunout,*)
-     s  nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys,
+     s  nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys,
      s  pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
      s  pctsrf(igout,is_sic)
@@ -3447 +3476 @@
       ENDIF
       
+!      first=.false.
 
       RETURN