source: LMDZ5/branches/LMDZ6_rc0/libf/phydev/physiq.F90

Last change on this file was 2160, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand

  • 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
12      USE dimphy, only : klon,klev
13      USE infotrac, only : nqtot
14      USE comgeomphy, only : rlatd
15      USE comcstphy, only : rg
16      USE iophy, only : histbeg_phy,histwrite_phy
17      USE ioipsl, only : getin,histvert,histdef,histend,ymds2ju
18      USE mod_phys_lmdz_para, only : jj_nb
19      USE phys_state_var_mod, only : phys_state_var_init
20
21#ifdef CPP_XIOS
22      USE xios, ONLY: xios_update_calendar
23      USE wxios, only: wxios_add_vaxis, wxios_set_timestep, wxios_closedef, &
24                       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
61integer,save :: itau=0 ! counter to count number of calls to physics
62!$OMP THREADPRIVATE(itau)
63real :: temp_newton(klon,klev)
64integer :: k
65logical, save :: first=.true.
66!$OMP THREADPRIVATE(first)
67
68! For I/Os
69integer :: itau0
70real :: zjulian
71real :: dtime
72integer :: nhori ! horizontal coordinate ID
73integer,save :: nid_hist ! output file ID
74!$OMP THREADPRIVATE(nid_hist)
75integer :: zvertid ! vertical coordinate ID
76integer,save :: iwrite_phys ! output every iwrite_phys physics step
77!$OMP THREADPRIVATE(iwrite_phys)
78integer,save :: iwrite_phys_omp ! intermediate variable to read iwrite_phys
79                                ! (must be shared by all threads)
80real :: t_ops ! frequency of the IOIPSL operations (eg average over...)
81real :: t_wrt ! frequency of the IOIPSL outputs
82
83! initializations
84if (debut) then ! Things to do only for the first call to physics
85! load initial conditions for physics (including the grid)
86  call phys_state_var_init() ! some initializations, required before calling phyetat0
87  call phyetat0("startphy.nc")
88
89! Initialize outputs:
90  itau0=0
91!$OMP MASTER
92  iwrite_phys_omp=1 !default: output every physics timestep
93  ! NB: getin() is not threadsafe; only one thread should call it.
94  call getin("iwrite_phys",iwrite_phys_omp)
95!$OMP END MASTER
96!$OMP BARRIER
97  iwrite_phys=iwrite_phys_omp
98  t_ops=pdtphys*iwrite_phys ! frequency of the IOIPSL operation
99  t_wrt=pdtphys*iwrite_phys ! frequency of the outputs in the file
100  ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0
101  !CALL ymds2ju(annee0, month, dayref, hour, zjulian)
102  call ymds2ju(1979, 1, 1, 0.0, zjulian)
103  dtime=pdtphys
104#ifndef CPP_IOIPSL_NO_OUTPUT
105  ! Initialize IOIPSL output file
106  call histbeg_phy("histins.nc",itau0,zjulian,dtime,nhori,nid_hist)
107#endif
108
109!$OMP MASTER
110
111#ifndef CPP_IOIPSL_NO_OUTPUT
112! IOIPSL
113  ! define vertical coordinate
114  call histvert(nid_hist,"presnivs","Vertical levels","Pa",klev, &
115                presnivs,zvertid,'down')
116  ! define variables which will be written in "histins.nc" file
117  call histdef(nid_hist,'temperature','Atmospheric temperature','K', &
118               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
119               'inst(X)',t_ops,t_wrt)
120  call histdef(nid_hist,'u','Eastward Zonal Wind','m/s', &
121               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
122               'inst(X)',t_ops,t_wrt)
123  call histdef(nid_hist,'v','Northward Meridional Wind','m/s', &
124               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
125               'inst(X)',t_ops,t_wrt)
126  call histdef(nid_hist,'ps','Surface Pressure','Pa', &
127               iim,jj_nb,nhori,1,1,1,zvertid,32, &
128               'inst(X)',t_ops,t_wrt)
129  ! end definition sequence
130  call histend(nid_hist)
131#endif
132
133#ifdef CPP_XIOS
134!XIOS
135    ! Declare available vertical axes to be used in output files:   
136    !CALL wxios_add_vaxis("presnivs", "dummy-not-used", klev, presnivs)
137    CALL wxios_add_vaxis("presnivs", klev, presnivs)
138
139    ! Declare time step length (in s):
140    CALL wxios_set_timestep(dtime)
141
142    !Finalize the context:
143    CALL wxios_closedef()
144#endif
145!$OMP END MASTER
146endif ! of if (debut)
147
148! increment local time counter itau
149itau=itau+1
150
151! set all tendencies to zero
152d_u(1:klon,1:klev)=0.
153d_v(1:klon,1:klev)=0.
154d_t(1:klon,1:klev)=0.
155d_qx(1:klon,1:klev,1:nqtot)=0.
156d_ps(1:klon)=0.
157
158! compute tendencies to return to the dynamics:
159! "friction" on the first layer
160d_u(1:klon,1)=-u(1:klon,1)/86400.
161d_v(1:klon,1)=-v(1:klon,1)/86400.
162! newtonian relaxation towards temp_newton()
163do k=1,klev
164  temp_newton(1:klon,k)=280.+cos(rlatd(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3
165  d_t(1:klon,k)=(temp_newton(1:klon,k)-t(1:klon,k))/1.e5
166enddo
167
168
169print*,'PHYDEV: itau=',itau
170
171! write some outputs:
172! IOIPSL
173#ifndef CPP_IOIPSL_NO_OUTPUT
174if (modulo(itau,iwrite_phys)==0) then
175  call histwrite_phy(nid_hist,.false.,"temperature",itau,t)
176  call histwrite_phy(nid_hist,.false.,"u",itau,u)
177  call histwrite_phy(nid_hist,.false.,"v",itau,v)
178  call histwrite_phy(nid_hist,.false.,"ps",itau,paprs(:,1))
179endif
180#endif
181
182!XIOS
183#ifdef CPP_XIOS
184!$OMP MASTER
185    !Increment XIOS time
186    CALL xios_update_calendar(itau)
187!$OMP END MASTER
188!$OMP BARRIER
189
190    !Send fields to XIOS: (NB these fields must also be defined as
191    ! <field id="..." /> in iodef.xml to be correctly used
192    CALL histwrite_phy("temperature",t)
193    CALL histwrite_phy("temp_newton",temp_newton)
194    CALL histwrite_phy("u",u)
195    CALL histwrite_phy("v",v)
196    CALL histwrite_phy("ps",paprs(:,1))
197#endif
198
199! if lastcall, then it is time to write "restartphy.nc" file
200if (lafin) then
201  call phyredem("restartphy.nc")
202endif
203
204end subroutine physiq
Note: See TracBrowser for help on using the repository browser.