source: LMDZ5/branches/LMDZ5-DOFOCO/libf/phydev/physiq.F90 @ 3893

Last change on this file since 3893 was 1759, checked in by Ehouarn Millour, 11 years ago

IOIPSL routine getin is not threadsafe, so when running in OpenMP, it should be called by only one thread (and the result copied to other threads in the case of threadprivate variables).
EM

File size: 6.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      IMPLICIT none
23#include "dimensions.h"
24
25      integer,parameter :: jjmp1=jjm+1-1/jjm
26      integer,parameter :: iip1=iim+1
27!
28! Routine argument:
29!
30      integer,intent(in) :: nlon ! number of atmospheric colums
31      integer,intent(in) :: nlev ! number of vertical levels (should be =klev)
32      real,intent(in) :: jD_cur ! current day number (Julian day)
33      real,intent(in) :: jH_cur ! current time of day (as fraction of day)
34      logical,intent(in) :: debut ! signals first call to physics
35      logical,intent(in) :: lafin ! signals last call to physics
36      real,intent(in) :: pdtphys ! physics time step (s)
37      real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa)
38      real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa)
39      real,intent(in) :: pphi(klon,klev) ! geopotential at mid-layer
40      real,intent(in) :: pphis(klon) ! surface geopotential
41      real,intent(in) :: presnivs(klev) ! pseudo-pressure (Pa) of mid-layers
42      integer,parameter :: longcles=20
43      real,intent(in) :: clesphy0(longcles) ! Not used
44      real,intent(in) :: u(klon,klev) ! eastward zonal wind (m/s)
45      real,intent(in) :: v(klon,klev) ! northward meridional wind (m/s)
46      real,intent(in) :: t(klon,klev) ! temperature (K)
47      real,intent(in) :: qx(klon,klev,nqtot) ! tracers (.../kg_air)
48      real,intent(in) :: flxmass_w(klon,klev) ! vertical mass flux
49      real,intent(out) :: d_u(klon,klev) ! physics tendency on u (m/s/s)
50      real,intent(out) :: d_v(klon,klev) ! physics tendency on v (m/s/s)
51      real,intent(out) :: d_t(klon,klev) ! physics tendency on t (K/s)
52      real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers
53      real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure
54      real,intent(in) :: dudyn(iim+1,jjmp1,klev) ! Not used
55!FH! REAL PVteta(klon,nbteta)
56!      REAL PVteta(klon,1)
57      real,intent(in) :: PVteta(klon,3) ! Not used ; should match definition
58                                        ! in calfis.F
59
60integer,save :: itau=0 ! counter to count number of calls to physics
61!$OMP THREADPRIVATE(itau)
62real :: temp_newton(klon,klev)
63integer :: k
64logical, save :: first=.true.
65!$OMP THREADPRIVATE(first)
66
67! For I/Os
68integer :: itau0
69real :: zjulian
70real :: dtime
71integer :: nhori ! horizontal coordinate ID
72integer,save :: nid_hist ! output file ID
73!$OMP THREADPRIVATE(nid_hist)
74integer :: zvertid ! vertical coordinate ID
75integer,save :: iwrite_phys ! output every iwrite_phys physics step
76!$OMP THREADPRIVATE(iwrite_phys)
77integer :: iwrite_phys_omp ! intermediate variable to read iwrite_phys
78real :: t_ops ! frequency of the IOIPSL operations (eg average over...)
79real :: t_wrt ! frequency of the IOIPSL outputs
80
81! initializations
82if (debut) then ! Things to do only for the first call to physics
83! load initial conditions for physics (including the grid)
84  call phys_state_var_init() ! some initializations, required before calling phyetat0
85  call phyetat0("startphy.nc")
86
87! Initialize outputs:
88  itau0=0
89!$OMP MASTER
90  iwrite_phys_omp=1 !default: output every physics timestep
91  ! NB: getin() is not threadsafe; only one thread should call it.
92  call getin("iwrite_phys",iwrite_phys_omp)
93!$OMP END MASTER
94!$OMP BARRIER
95  iwrite_phys=iwrite_phys_omp
96  t_ops=pdtphys*iwrite_phys ! frequency of the IOIPSL operation
97  t_wrt=pdtphys*iwrite_phys ! frequency of the outputs in the file
98  ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0
99  !CALL ymds2ju(annee0, month, dayref, hour, zjulian)
100  call ymds2ju(1979, 1, 1, 0.0, zjulian)
101  dtime=pdtphys
102  call histbeg_phy("histins.nc",itau0,zjulian,dtime,nhori,nid_hist)
103!$OMP MASTER
104  ! define vertical coordinate
105  call histvert(nid_hist,"presnivs","Vertical levels","Pa",klev, &
106                presnivs,zvertid,'down')
107  ! define variables which will be written in "histins.nc" file
108  call histdef(nid_hist,'temperature','Atmospheric temperature','K', &
109               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
110               'inst(X)',t_ops,t_wrt)
111  call histdef(nid_hist,'u','Eastward Zonal Wind','m/s', &
112               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
113               'inst(X)',t_ops,t_wrt)
114  call histdef(nid_hist,'v','Northward Meridional Wind','m/s', &
115               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
116               'inst(X)',t_ops,t_wrt)
117  call histdef(nid_hist,'ps','Surface Pressure','Pa', &
118               iim,jj_nb,nhori,1,1,1,zvertid,32, &
119               'inst(X)',t_ops,t_wrt)
120  ! end definition sequence
121  call histend(nid_hist)
122!$OMP END MASTER
123endif ! of if (debut)
124
125! increment counter itau
126itau=itau+1
127
128! set all tendencies to zero
129d_u(1:klon,1:klev)=0.
130d_v(1:klon,1:klev)=0.
131d_t(1:klon,1:klev)=0.
132d_qx(1:klon,1:klev,1:nqtot)=0.
133d_ps(1:klon)=0.
134
135! compute tendencies to return to the dynamics:
136! "friction" on the first layer
137d_u(1:klon,1)=-u(1:klon,1)/86400.
138d_v(1:klon,1)=-v(1:klon,1)/86400.
139! newtonian relaxation towards temp_newton()
140do k=1,klev
141  temp_newton(1:klon,k)=280.+cos(rlatd(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3
142  d_t(1:klon,k)=(temp_newton(1:klon,k)-t(1:klon,k))/1.e5
143enddo
144
145
146!print*,'PHYDEV: itau=',itau
147
148! write some outputs:
149if (modulo(itau,iwrite_phys)==0) then
150  call histwrite_phy(nid_hist,.false.,"temperature",itau,t)
151  call histwrite_phy(nid_hist,.false.,"u",itau,u)
152  call histwrite_phy(nid_hist,.false.,"v",itau,v)
153  call histwrite_phy(nid_hist,.false.,"ps",itau,paprs(:,1))
154endif
155
156! if lastcall, then it is time to write "restartphy.nc" file
157if (lafin) then
158  call phyredem("restartphy.nc")
159endif
160
161end
Note: See TracBrowser for help on using the repository browser.