source: LMDZ5/trunk/libf/phydev/physiq.F90 @ 1919

Last change on this file since 1919 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 7.1 KB
Line 
1! $Id: physiq.F 1565 2011-08-31 12:53:29Z jghattas $
2!#define IO_DEBUG
3
4      SUBROUTINE physiq (nlon,nlev, &
5     &            debut,lafin,jD_cur, jH_cur,pdtphys, &
6     &            paprs,pplay,pphi,pphis,presnivs,clesphy0, &
7     &            u,v,t,qx, &
8     &            flxmass_w, &
9     &            d_u, d_v, d_t, d_qx, d_ps &
10     &            , dudyn &
11     &            , PVteta)
12
13      USE dimphy, only : klon,klev
14      USE infotrac, only : nqtot
15      USE comgeomphy, only : rlatd
16      USE comcstphy, only : rg
17      USE iophy, only : histbeg_phy,histwrite_phy
18      USE ioipsl, only : getin,histvert,histdef,histend,ymds2ju
19      USE mod_phys_lmdz_para, only : jj_nb
20      USE phys_state_var_mod, only : phys_state_var_init
21
22#ifdef CPP_XIOS
23      USE wxios, only: wxios_add_vaxis, wxios_set_timestep, wxios_closedef, &
24                       wxios_update_calendar, histwrite_phy
25#endif
26
27      IMPLICIT none
28#include "dimensions.h"
29
30      integer,parameter :: jjmp1=jjm+1-1/jjm
31      integer,parameter :: iip1=iim+1
32!
33! Routine argument:
34!
35      integer,intent(in) :: nlon ! number of atmospheric colums
36      integer,intent(in) :: nlev ! number of vertical levels (should be =klev)
37      real,intent(in) :: jD_cur ! current day number (Julian day)
38      real,intent(in) :: jH_cur ! current time of day (as fraction of day)
39      logical,intent(in) :: debut ! signals first call to physics
40      logical,intent(in) :: lafin ! signals last call to physics
41      real,intent(in) :: pdtphys ! physics time step (s)
42      real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa)
43      real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa)
44      real,intent(in) :: pphi(klon,klev) ! geopotential at mid-layer
45      real,intent(in) :: pphis(klon) ! surface geopotential
46      real,intent(in) :: presnivs(klev) ! pseudo-pressure (Pa) of mid-layers
47      integer,parameter :: longcles=20
48      real,intent(in) :: clesphy0(longcles) ! Not used
49      real,intent(in) :: u(klon,klev) ! eastward zonal wind (m/s)
50      real,intent(in) :: v(klon,klev) ! northward meridional wind (m/s)
51      real,intent(in) :: t(klon,klev) ! temperature (K)
52      real,intent(in) :: qx(klon,klev,nqtot) ! tracers (.../kg_air)
53      real,intent(in) :: flxmass_w(klon,klev) ! vertical mass flux
54      real,intent(out) :: d_u(klon,klev) ! physics tendency on u (m/s/s)
55      real,intent(out) :: d_v(klon,klev) ! physics tendency on v (m/s/s)
56      real,intent(out) :: d_t(klon,klev) ! physics tendency on t (K/s)
57      real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers
58      real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure
59      real,intent(in) :: dudyn(iim+1,jjmp1,klev) ! Not used
60!FH! REAL PVteta(klon,nbteta)
61!      REAL PVteta(klon,1)
62      real,intent(in) :: PVteta(klon,3) ! Not used ; should match definition
63                                        ! in calfis.F
64
65integer,save :: itau=0 ! counter to count number of calls to physics
66!$OMP THREADPRIVATE(itau)
67real :: temp_newton(klon,klev)
68integer :: k
69logical, save :: first=.true.
70!$OMP THREADPRIVATE(first)
71
72! For I/Os
73integer :: itau0
74real :: zjulian
75real :: dtime
76integer :: nhori ! horizontal coordinate ID
77integer,save :: nid_hist ! output file ID
78!$OMP THREADPRIVATE(nid_hist)
79integer :: zvertid ! vertical coordinate ID
80integer,save :: iwrite_phys ! output every iwrite_phys physics step
81!$OMP THREADPRIVATE(iwrite_phys)
82integer,save :: iwrite_phys_omp ! intermediate variable to read iwrite_phys
83                                ! (must be shared by all threads)
84real :: t_ops ! frequency of the IOIPSL operations (eg average over...)
85real :: t_wrt ! frequency of the IOIPSL outputs
86
87! initializations
88if (debut) then ! Things to do only for the first call to physics
89! load initial conditions for physics (including the grid)
90  call phys_state_var_init() ! some initializations, required before calling phyetat0
91  call phyetat0("startphy.nc")
92
93! Initialize outputs:
94  itau0=0
95!$OMP MASTER
96  iwrite_phys_omp=1 !default: output every physics timestep
97  ! NB: getin() is not threadsafe; only one thread should call it.
98  call getin("iwrite_phys",iwrite_phys_omp)
99!$OMP END MASTER
100!$OMP BARRIER
101  iwrite_phys=iwrite_phys_omp
102  t_ops=pdtphys*iwrite_phys ! frequency of the IOIPSL operation
103  t_wrt=pdtphys*iwrite_phys ! frequency of the outputs in the file
104  ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0
105  !CALL ymds2ju(annee0, month, dayref, hour, zjulian)
106  call ymds2ju(1979, 1, 1, 0.0, zjulian)
107  dtime=pdtphys
108#ifndef CPP_NO_IOIPSL
109  ! Initialize IOIPSL output file
110  call histbeg_phy("histins.nc",itau0,zjulian,dtime,nhori,nid_hist)
111#endif
112
113!$OMP MASTER
114
115#ifndef CPP_NO_IOIPSL
116! IOIPSL
117  ! define vertical coordinate
118  call histvert(nid_hist,"presnivs","Vertical levels","Pa",klev, &
119                presnivs,zvertid,'down')
120  ! define variables which will be written in "histins.nc" file
121  call histdef(nid_hist,'temperature','Atmospheric temperature','K', &
122               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
123               'inst(X)',t_ops,t_wrt)
124  call histdef(nid_hist,'u','Eastward Zonal Wind','m/s', &
125               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
126               'inst(X)',t_ops,t_wrt)
127  call histdef(nid_hist,'v','Northward Meridional Wind','m/s', &
128               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
129               'inst(X)',t_ops,t_wrt)
130  call histdef(nid_hist,'ps','Surface Pressure','Pa', &
131               iim,jj_nb,nhori,1,1,1,zvertid,32, &
132               'inst(X)',t_ops,t_wrt)
133  ! end definition sequence
134  call histend(nid_hist)
135#endif
136
137#ifdef CPP_XIOS
138!XIOS
139    ! Déclaration de l'axe vertical du fichier:   
140    CALL wxios_add_vaxis("presnivs", "histins", klev, presnivs)
141
142    !Déclaration du pas de temps:
143    CALL wxios_set_timestep(dtime)
144
145    !Finalisation du contexte:
146    CALL wxios_closedef()
147#endif
148!$OMP END MASTER
149endif ! of if (debut)
150
151! increment local time counter itau
152itau=itau+1
153
154! set all tendencies to zero
155d_u(1:klon,1:klev)=0.
156d_v(1:klon,1:klev)=0.
157d_t(1:klon,1:klev)=0.
158d_qx(1:klon,1:klev,1:nqtot)=0.
159d_ps(1:klon)=0.
160
161! compute tendencies to return to the dynamics:
162! "friction" on the first layer
163d_u(1:klon,1)=-u(1:klon,1)/86400.
164d_v(1:klon,1)=-v(1:klon,1)/86400.
165! newtonian relaxation towards temp_newton()
166do k=1,klev
167  temp_newton(1:klon,k)=280.+cos(rlatd(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3
168  d_t(1:klon,k)=(temp_newton(1:klon,k)-t(1:klon,k))/1.e5
169enddo
170
171
172print*,'PHYDEV: itau=',itau
173
174! write some outputs:
175! IOIPSL
176#ifndef CPP_NO_IOIPSL
177if (modulo(itau,iwrite_phys)==0) then
178  call histwrite_phy(nid_hist,.false.,"temperature",itau,t)
179  call histwrite_phy(nid_hist,.false.,"u",itau,u)
180  call histwrite_phy(nid_hist,.false.,"v",itau,v)
181  call histwrite_phy(nid_hist,.false.,"ps",itau,paprs(:,1))
182endif
183#endif
184
185!XIOS
186#ifdef CPP_XIOS
187!$OMP MASTER
188    !Increment XIOS time
189    CALL wxios_update_calendar(itau)
190!$OMP END MASTER
191!$OMP BARRIER
192
193    !Send fields to XIOS:
194    CALL histwrite_phy("temperature",t)
195    CALL histwrite_phy("u",u)
196    CALL histwrite_phy("v",v)
197    CALL histwrite_phy("ps",paprs(:,1))
198#endif
199
200! if lastcall, then it is time to write "restartphy.nc" file
201if (lafin) then
202  call phyredem("restartphy.nc")
203endif
204
205end subroutine physiq
Note: See TracBrowser for help on using the repository browser.