Changeset 1707 for LMDZ5/branches/testing/libf/phydev/physiq.F90

Timestamp:

Jan 11, 2013, 10:19:19 AM (11 years ago)

Author:

Laurent Fairhead

Message:

Version testing basée sur la r1706

---

Testing release based on r1706

Location:

LMDZ5/branches/testing

Files:

• . (modified) (1 prop)
• libf/phydev/physiq.F90 (modified) (1 diff)

LMDZ5/branches/testing/libf/phydev/physiq.F90

@@ -11 +11 @@
      &            , PVteta)
 
-      USE dimphy
-      USE infotrac
-      USE comgeomphy
+      USE dimphy, only : klon,klev
+      USE infotrac, only : nqtot
+      USE comgeomphy, only : rlatd
+      USE comcstphy, only : rg
+      USE iophy, only : histbeg_phy,histwrite_phy
+      USE ioipsl, only : getin,histvert,histdef,histend,ymds2ju
+      USE mod_phys_lmdz_para, only : jj_nb
+      USE phys_state_var_mod, only : phys_state_var_init
 
       IMPLICIT none
-!======================================================================
-! Objet: Moniteur general de la physique du modele
-!======================================================================
+#include "dimensions.h"
+
+      integer,parameter :: jjmp1=jjm+1-1/jjm
+      integer,parameter :: iip1=iim+1
 !
-!  Arguments:
+! Routine argument:
 !
-! nlon----input-I-nombre de points horizontaux
-! nlev----input-I-nombre de couches verticales, doit etre egale a klev
-! debut---input-L-variable logique indiquant le premier passage
-! lafin---input-L-variable logique indiquant le dernier passage
-! jD_cur       -R-jour courant a l'appel de la physique (jour julien)
-! jH_cur       -R-heure courante a l'appel de la physique (jour julien)
-! pdtphys-input-R-pas d'integration pour la physique (seconde)
-! paprs---input-R-pression pour chaque inter-couche (en Pa)
-! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
-! pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
-! pphis---input-R-geopotentiel du sol
-! presnivs-input_R_pressions approximat. des milieux couches ( en PA)
-! u-------input-R-vitesse dans la direction X (de O a E) en m/s
-! v-------input-R-vitesse Y (de S a N) en m/s
-! t-------input-R-temperature (K)
-! qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
-! d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
-! d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
-! flxmass_w -input-R- flux de masse verticale
-! d_u-----output-R-tendance physique de "u" (m/s/s)
-! d_v-----output-R-tendance physique de "v" (m/s/s)
-! d_t-----output-R-tendance physique de "t" (K/s)
-! d_qx----output-R-tendance physique de "qx" (kg/kg/s)
-! d_ps----output-R-tendance physique de la pression au sol
-!IM
-! PVteta--output-R-vorticite potentielle a des thetas constantes
-!======================================================================
-#include "dimensions.h"
-#include "comcstphy.h"
+      integer,intent(in) :: nlon ! number of atmospheric colums
+      integer,intent(in) :: nlev ! number of vertical levels (should be =klev)
+      real,intent(in) :: jD_cur ! current day number (Julian day)
+      real,intent(in) :: jH_cur ! current time of day (as fraction of day)
+      logical,intent(in) :: debut ! signals first call to physics
+      logical,intent(in) :: lafin ! signals last call to physics
+      real,intent(in) :: pdtphys ! physics time step (s)
+      real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa)
+      real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa)
+      real,intent(in) :: pphi(klon,klev) ! geopotential at mid-layer
+      real,intent(in) :: pphis(klon) ! surface geopotential
+      real,intent(in) :: presnivs(klev) ! pseudo-pressure (Pa) of mid-layers
+      integer,parameter :: longcles=20
+      real,intent(in) :: clesphy0(longcles) ! Not used
+      real,intent(in) :: u(klon,klev) ! eastward zonal wind (m/s)
+      real,intent(in) :: v(klon,klev) ! northward meridional wind (m/s)
+      real,intent(in) :: t(klon,klev) ! temperature (K)
+      real,intent(in) :: qx(klon,klev,nqtot) ! tracers (.../kg_air)
+      real,intent(in) :: flxmass_w(klon,klev) ! vertical mass flux
+      real,intent(out) :: d_u(klon,klev) ! physics tendency on u (m/s/s)
+      real,intent(out) :: d_v(klon,klev) ! physics tendency on v (m/s/s)
+      real,intent(out) :: d_t(klon,klev) ! physics tendency on t (K/s)
+      real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers
+      real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure
+      real,intent(in) :: dudyn(iim+1,jjmp1,klev) ! Not used
+!FH! REAL PVteta(klon,nbteta)
+!      REAL PVteta(klon,1)
+      real,intent(in) :: PVteta(klon,3) ! Not used ; should match definition
+                                        ! in calfis.F
 
-      integer jjmp1
-      parameter (jjmp1=jjm+1-1/jjm)
-      integer iip1
-      parameter (iip1=iim+1)
+integer,save :: itau=0 ! counter to count number of calls to physics
+!$OMP THREADPRIVATE(itau)
+real :: temp_newton(klon,klev)
+integer :: k
+logical, save :: first=.true.
+!$OMP THREADPRIVATE(first)
 
-      INTEGER ivap          ! indice de traceurs pour vapeur d'eau
-      PARAMETER (ivap=1)
-      INTEGER iliq          ! indice de traceurs pour eau liquide
-      PARAMETER (iliq=2)
+! For I/Os
+integer :: itau0
+real :: zjulian
+real :: dtime
+integer :: nhori ! horizontal coordinate ID
+integer,save :: nid_hist ! output file ID
+!$OMP THREADPRIVATE(nid_hist)
+integer :: zvertid ! vertical coordinate ID
+integer,save :: iwrite_phys ! output every iwrite_phys physics step
+!$OMP THREADPRIVATE(iwrite_phys)
+real :: t_ops ! frequency of the IOIPSL operations (eg average over...)
+real :: t_wrt ! frequency of the IOIPSL outputs
 
-!
-!
-! Variables argument:
-!
-      INTEGER nlon
-      INTEGER nlev
-      REAL, intent(in):: jD_cur, jH_cur
+! initializations
+if (debut) then ! Things to do only for the first call to physics 
+! load initial conditions for physics (including the grid)
+  call phys_state_var_init() ! some initializations, required before calling phyetat0
+  call phyetat0("startphy.nc")
 
-      REAL pdtphys
-      LOGICAL debut, lafin
-      REAL paprs(klon,klev+1)
-      REAL pplay(klon,klev)
-      REAL pphi(klon,klev)
-      REAL pphis(klon)
-      REAL presnivs(klev)
-      REAL znivsig(klev)
-      real pir
+! Initialize outputs:
+  itau0=0
+  iwrite_phys=1 !default: output every physics timestep
+  call getin("iwrite_phys",iwrite_phys)
+  t_ops=pdtphys*iwrite_phys ! frequency of the IOIPSL operation
+  t_wrt=pdtphys*iwrite_phys ! frequency of the outputs in the file
+  ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0
+  !CALL ymds2ju(annee0, month, dayref, hour, zjulian)
+  call ymds2ju(1979, 1, 1, 0.0, zjulian)
+  dtime=pdtphys
+  call histbeg_phy("histins.nc",itau0,zjulian,dtime,nhori,nid_hist)
+!$OMP MASTER
+  ! define vertical coordinate
+  call histvert(nid_hist,"presnivs","Vertical levels","Pa",klev, &
+                presnivs,zvertid,'down')
+  ! define variables which will be written in "histins.nc" file
+  call histdef(nid_hist,'temperature','Atmospheric temperature','K', &
+               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
+               'inst(X)',t_ops,t_wrt)
+  call histdef(nid_hist,'u','Eastward Zonal Wind','m/s', &
+               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
+               'inst(X)',t_ops,t_wrt)
+  call histdef(nid_hist,'v','Northward Meridional Wind','m/s', &
+               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
+               'inst(X)',t_ops,t_wrt)
+  call histdef(nid_hist,'ps','Surface Pressure','Pa', &
+               iim,jj_nb,nhori,1,1,1,zvertid,32, &
+               'inst(X)',t_ops,t_wrt)
+  ! end definition sequence
+  call histend(nid_hist)
+!$OMP END MASTER
+endif ! of if (debut)
 
-      REAL u(klon,klev)
-      REAL v(klon,klev)
-      REAL t(klon,klev),theta(klon,klev)
-      REAL qx(klon,klev,nqtot)
-      REAL flxmass_w(klon,klev)
-      REAL omega(klon,klev) ! vitesse verticale en Pa/s
-      REAL d_u(klon,klev)
-      REAL d_v(klon,klev)
-      REAL d_t(klon,klev)
-      REAL d_qx(klon,klev,nqtot)
-      REAL d_ps(klon)
-      real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
-!IM definition dynamique o_trac dans phys_output_open
-!      type(ctrl_out) :: o_trac(nqtot)
-!FH! REAL PVteta(klon,nbteta)
-      REAL PVteta(klon,1)
-      REAL dudyn(iim+1,jjmp1,klev)
+! increment counter itau
+itau=itau+1
 
-    INTEGER        longcles
-    PARAMETER    ( longcles = 20 )
+! set all tendencies to zero
+d_u(1:klon,1:klev)=0.
+d_v(1:klon,1:klev)=0.
+d_t(1:klon,1:klev)=0.
+d_qx(1:klon,1:klev,1:nqtot)=0.
+d_ps(1:klon)=0.
 
-real temp_newton(klon,klev)
-integer k
-logical, save :: first=.true.
-
-      REAL clesphy0( longcles      )
-
-d_u=0.
-d_v=0.
-d_t=0.
-d_qx=0.
-d_ps=0.
-
-     d_u(:,1)=-u(:,1)/86400.
-     do k=1,klev
-        temp_newton(:,k)=280.+cos(rlatd(:))*40.-pphi(:,k)/rg*6.e-3
-        d_t(:,k)=(temp_newton(:,k)-t(:,k))/1.e5
-     enddo
+! compute tendencies to return to the dynamics:
+! "friction" on the first layer
+d_u(1:klon,1)=-u(1:klon,1)/86400.
+d_v(1:klon,1)=-v(1:klon,1)/86400.
+! newtonian relaxation towards temp_newton()
+do k=1,klev
+  temp_newton(1:klon,k)=280.+cos(rlatd(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3
+  d_t(1:klon,k)=(temp_newton(1:klon,k)-t(1:klon,k))/1.e5
+enddo
 
 
-      print*,'COUCOU PHYDEV'
-      return
-      end
+!print*,'PHYDEV: itau=',itau
+
+! write some outputs:
+if (modulo(itau,iwrite_phys)==0) then
+  call histwrite_phy(nid_hist,.false.,"temperature",itau,t)
+  call histwrite_phy(nid_hist,.false.,"u",itau,u)
+  call histwrite_phy(nid_hist,.false.,"v",itau,v)
+  call histwrite_phy(nid_hist,.false.,"ps",itau,paprs(:,1))
+endif
+
+! if lastcall, then it is time to write "restartphy.nc" file
+if (lafin) then
+  call phyredem("restartphy.nc")
+endif
+
+end