source: LMDZ6/branches/Portage_acc/libf/dyn3d/driver.F90 @ 4166

Last change on this file since 4166 was 4166, checked in by Ehouarn Millour, 2 years ago

Update driver to test repeated calling of physics.
Add some intent() to ajsec.
And !$acc contructs for physiq arguments.
EM

File size: 15.6 KB
Line 
1PROGRAM driver
2
3  USE IOIPSL, ONLY: getin, ymds2ju, ioconf_calendar
4  USE mod_const_mpi, ONLY: COMM_LMDZ
5  USE control_mod, ONLY: planet_type, config_inca, offline
6  USE control_mod, ONLY: nday, day_step, iphysiq, dayref, anneeref
7  USE infotrac, ONLY: type_trac, nqtot, infotrac_init
8  USE temps_mod, ONLY: calend,start_time,annee_ref,day_ref
9  USE temps_mod, ONLY: itau_dyn,itau_phy,day_ini,jD_ref,jH_ref,day_end,year_len
10  USE comconst_mod, ONLY: cpp, kappa, daysec, dtphys, dtvr, g, r, rad
11  USE comconst_mod, ONLY: daylen, year_day, pi, omeg
12  USE logic_mod, ONLY: iflag_phys
13  USE comvert_mod, ONLY: preff, pa
14  USE inigeomphy_mod, ONLY: inigeomphy
15  USE iniphysiq_mod, ONLY: iniphysiq
16
17  IMPLICIT NONE
18
19  include "dimensions.h"
20  include "paramet.h"
21  include "comgeom.h"
22
23!  INTEGER :: step_per_day = 24 ! # of physics steps per day
24
25  INTEGER :: ngrid ! # of grid points on physics grid = 2+(jjm-1)*iim-1/jjm
26  INTEGER :: ig,l
27
28!  REAL, PARAMETER :: unjours=86400.,    & ! solar day in seconds
29!       &             radius=6.4e6,      & ! planetary radius
30!       &             g=9.8,             & ! gravity
31!       &             cpp=1004.,         & ! Cp
32!       &             kappa=296.945007/cpp ! R/Cp
33!  REAL :: timestep !=unjours/step_per_day ! physics time step (s)
34
35!  REAL :: psurf=1e5, ptop=1e4, Temp=250. ! initial values of surface pressure, temperature
36
37  ! dynamics (to load data from start.nc)
38  REAL,ALLOCATABLE :: ucov(:,:) ! covariant zonal wind
39  REAL,ALLOCATABLE :: vcov(:,:) ! covariant meridional wind
40  REAL,ALLOCATABLE :: teta(:,:) ! potential temperature
41  REAL,ALLOCATABLE :: q(:,:,:) ! tracers
42  REAL,ALLOCATABLE :: masse(:,:) ! air mass
43  REAL,ALLOCATABLE :: ps(:) ! surface pressure
44  REAL,ALLOCATABLE :: phis(:) ! surface geopotential
45  REAL :: time_0
46 
47  ! physics
48  REAL,ALLOCATABLE :: pplev(:,:)  ! pressure at interfaces
49  REAL,ALLOCATABLE :: pplay(:,:) ! pressure at full levels
50  REAL,ALLOCATABLE :: pphi(:,:) ! geopotential at full levels
51  REAL,ALLOCATABLE :: pphis(:) ! surface geopotential
52  REAL,ALLOCATABLE :: pt(:,:) ! temperature at full levels
53  REAL,ALLOCATABLE :: pu(:,:) ! zonal wind at full levels
54  REAL,ALLOCATABLE :: pv(:,:) ! meridional wind at full levels
55  REAL,ALLOCATABLE :: pq(:,:,:)! tracers
56  REAL,ALLOCATABLE :: pr(:,:) ! vorticity
57
58
59! 0. Preliminary stuff: read in parameters from run.def
60  call conf_gcm( 99, .TRUE.)
61!  call conf_planete
62 
63  daysec=86400. ! day length(s)
64  daylen=daysec
65  year_day=360 ! # days per year
66  rad=6.4e6 ! planet radius (m)
67  g=9.81 ! gravity
68  cpp=1004. ! Cp
69  r=287.05967 ! R
70  kappa=287.05967/cpp ! R/Cp
71  preff=101325.0
72  pa=preff/2.
73  pi=2.*asin(1.)
74  omeg=2.*pi/daysec*(1./daylen+1./year_day)
75 
76!  nday = 1 ! simulation lenght (days)
77!  CALL getin("nday",nday)
78  WRITE(*,*)"driver: nday=",nday
79 
80!  CALL getin("day_step",day_step)
81  WRITE(*,*)"driver: day_step=",day_step
82  dtphys=daysec/(day_step/iphysiq) ! physics time step (s)
83  dtvr=daysec/day_step
84  WRITE(*,*)"driver: timestep= dtphys=",dtphys
85 
86 
87 
88! 1. Initialize setup (geometry, physics)
89! 1.0. Parameters
90  planet_type="earth"
91  COMM_LMDZ=1
92  iflag_phys=1
93 
94
95! 1.1. Tracers-related
96  type_trac="lmdz"
97  config_inca="none"
98  offline=.false.
99 
100  call infotrac_init
101 
102  ! 1.2. Geometry
103
104!  CALL init_latlon(lat, lon)
105!  WRITE(*,*)"driver: LAT = " , minval(lat), maxval(lat)
106!  WRITE(*,*)"driver: LON = " , minval(lon), maxval(lon)
107  ! initialize dyn constants and geometry
108  CALL iniconst
109  CALL inigeom
110
111  ! 1.3. Initial conditions
112  ! dynamics
113  ALLOCATE(ucov(ip1jmp1,llm))
114  ALLOCATE(vcov(ip1jm,llm))
115  ALLOCATE(teta(ip1jmp1,llm))
116  ALLOCATE(q(ip1jmp1,llm,nqtot))
117  ALLOCATE(masse(ip1jmp1,llm))
118  ALLOCATE(ps(ip1jmp1))
119  ALLOCATE(phis(ip1jmp1))
120 
121  CALL dynetat0("start.nc",vcov,ucov,teta,q,masse,ps,phis,time_0)
122 
123
124  ! 1.3. Calendar
125!  annee_ref=2000
126!  dayref=1
127  if (calend == 'earth_360d' ) then
128    call ioconf_calendar('360d')
129  else
130    write(*,*) "potential calendar problem with IOIPSL!!"
131  endif
132  year_len=360
133
134!  annee_ref = anneeref
135!  day_ref = dayref
136!  day_ini = dayref
137  itau_dyn = 0
138  itau_phy = 0
139
140  ! Calendar
141!  call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
142  call ymds2ju(annee_ref, 1, day_ref, 0., jD_ref)
143  jH_ref = jD_ref - int(jD_ref)
144  jD_ref = int(jD_ref)
145!  call ioconf_startdate(INT(jD_ref), jH_ref)
146 
147
148  day_end = day_ini + nday
149
150  ! 1.4. Initialize physics
151  ngrid=2+(jjm-1)*iim-1/jjm
152  write(*,*) "driver: ngrid=",ngrid
153  CALL iniphysiq(iim,jjm,llm, &
154          ngrid,comm_lmdz, &
155          daysec,day_ini,dtphys, &
156          rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
157          iflag_phys)
158
159 
160! 2. Initialize fields on physics grid
161
162  ALLOCATE(pplev(ngrid,llm+1)) ! pressure at interfaces
163  ALLOCATE(pplay(ngrid,llm)) ! pressure at full levels
164  ALLOCATE(pphi(ngrid,llm)) ! geopotential at full levels
165  ALLOCATE(pphis(ngrid)) ! surface geopotential
166  ALLOCATE(pt(ngrid,llm)) ! temperature (K) at full levels
167  ALLOCATE(pu(ngrid,llm)) ! zonal wind (m/s) at full levels
168  ALLOCATE(pr(ngrid,llm)) ! relative vorticity
169  ALLOCATE(pv(ngrid,llm)) ! meridional wind at full levels
170  ALLOCATE(pq(ngrid,llm,nqtot)) ! tracers
171
172  CALL init_temperature_pressure_geopot(ps,masse,teta,phis, &
173                                        ngrid,pplev,pplay,pphi,pphis,pt)
174  CALL init_winds_tracers(ucov,vcov,q,nqtot, &
175                          ngrid,pu,pv,pr,pq)
176
177! 3. Temporal loop 
178
179  ! First time on GPU
180  !$acc data create(pplev, pplay, pphi, pphis, pt, pu, pv, pq, pr)
181
182  CALL timeloop(ngrid,llm,nqtot,nday,day_step/iphysiq,dtphys,&
183                pplev,pplay,pphi,pphis,&
184                pu,pv,pr,pt,pq)
185
186  !$acc end data
187
188  ! output final state as plain text file
189  open(10,file="driver.dat")
190  do ig=1,ngrid
191   do l=1,llm
192    write(10,"(I5,I3,X,20(1PE15.8,X))") &
193         ig,l,pu(ig,l),pv(ig,l),pt(ig,l)!,pq(ig,l,1:nqtot)
194   enddo
195  enddo
196  close(10)
197
198  write(*,*)"driver: Everything is cool :-)"
199
200CONTAINS
201
202!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
203
204   SUBROUTINE timeloop(ngrid,llm,nqtot,nday,calls_per_day,dtphys,&
205                pplev,pplay,pphi,pphis,&
206                pu,pv,pr,pt,pq)
207
208   USE comvert_mod, ONLY: presnivs
209   USE physiq_mod, ONLY: physiq
210   USE comconst_mod, ONLY: r
211   IMPLICIT NONE
212
213   INTEGER,INTENT(IN) :: ngrid ! # number of atmospheric columns
214   INTEGER,INTENT(IN) :: llm ! # number of vertical levels
215   INTEGER,INTENT(IN) :: nqtot ! # of tracers
216   INTEGER,INTENT(IN) :: nday ! # of days to run
217   INTEGER,INTENT(IN) :: calls_per_day ! # number of calls to physics per day
218   REAL,INTENT(IN) :: dtphys ! physics time step (s)
219   REAL,INTENT(INOUT) :: pplev(ngrid,llm+1) ! pressure at interfaces
220   REAL,INTENT(INOUT) :: pplay(ngrid,llm) ! pressure at mid-layer
221   REAL,INTENT(INOUT) :: pphi(ngrid,llm) ! geopotential at mid-layer
222   REAL,INTENT(INOUT) :: pphis(ngrid) ! surface geopotential
223   REAL,INTENT(INOUT) :: pu(ngrid,llm) ! zonal wind (m/s)
224   REAL,INTENT(INOUT) :: pv(ngrid,llm) ! meridional wind (m/s)
225   REAL,INTENT(INOUT) :: pr(ngrid,llm) ! relative vorticity
226   REAL,INTENT(INOUT) :: pt(ngrid,llm) ! temperature (K)
227   REAL,INTENT(INOUT) :: pq(ngrid,llm,nqtot) ! tracers
228   
229   LOGICAL,SAVE :: firstcall=.true.
230   LOGICAL,SAVE :: lastcall=.false.
231   REAL :: JD_cur ! Julian day
232   REAL :: JH_cur ! Julian hour (fraction of day)
233   ! tendencies from the physics
234   REAL :: pdu(ngrid,llm)
235   REAL :: pdv(ngrid,llm)
236   REAL :: pdt(ngrid,llm)
237   REAL :: pdq(ngrid,llm,nqtot)
238   REAL :: pdps(ngrid)
239   
240   INTEGER :: iday,istep
241   INTEGER :: ig,l
242   REAL :: rho,gdz
243
244   REAL :: flxw(ngrid,llm) ! vertical mass flux ! set to zero here.
245   REAL :: phi_top(ngrid)
246   !$acc data present(pplev,pplay,pphi,pphis,pu,pv,pt,pq) &
247   !$acc & create(phi_top,flxw,pdu,pdv,pdt,pdq,pdps)
248   
249   flxw(1:ngrid,1:llm)=0.
250   
251   DO iday=1,nday
252     JD_cur=iday
253     DO istep = 0, calls_per_day-1
254       JH_cur=istep/calls_per_day
255       if ((iday==nday).and.(istep==calls_per_day-1)) lastcall=.true.
256       CALL physiq(ngrid,llm,firstcall,lastcall,dtphys, &
257                   pplev,pplay,pphi,pphis, &
258                   presnivs, &
259                   pu,pv,pr,pt,pq,flxw, &
260                   pdu,pdv,pdt,pdq,pdps)
261       firstcall=.false.
262       
263       ! add increments sent back from physics:
264       !$acc kernels default(present)
265       pt(:,:)=pt(:,:)+dtphys*pdt(:,:)
266       pu(:,:)=pu(:,:)+dtphys*pdu(:,:)
267       pv(:,:)=pv(:,:)+dtphys*pdv(:,:)
268       pq(:,:,:)=pq(:,:,:)+dtphys*pdq(:,:,:)
269       
270       ! Recomputation of relative geopotential:
271       phi_top(:)=0
272       !$acc loop private(rho, gdz)
273       do l=1,llm
274         do ig=1,ngrid
275           ! rho=p/RT
276           rho=r*pplay(ig,l)/pt(ig,l)
277           gdz=(pplev(ig,l)-pplev(ig,l+1))/rho ! layer thickness * g
278           pphi(ig,l)=phi_top(ig)+0.5*gdz ! geopotential at mid-layer
279           phi_top(ig)=phi_top(ig)+gdz ! geopotential at layer top
280         enddo
281       enddo
282       !$acc end kernels
283       
284     ENDDO ! of DO istep = 0, day_step-1
285   ENDDO ! of DO iday=0,nday-1
286   
287   !$acc end data
288   
289   END SUBROUTINE timeloop
290
291
292!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
293
294   SUBROUTINE init_winds_tracers(ucov,vcov,q,nqtot, &
295                                 ngrid,pu,pv,pr,pq)
296
297    USE comconst_mod, ONLY: pi
298   
299    IMPLICIT NONE
300
301    include "dimensions.h"
302    include "paramet.h"
303    include "comgeom2.h"
304
305     REAL,INTENT(IN) :: ucov(iip1,jjp1,llm) ! covariant meridional velocity
306     REAL,INTENT(IN) :: vcov(iip1,jjm,llm) ! covariant zonal wind
307     REAL,INTENT(IN) :: q(iip1,jjp1,llm,nqtot) ! tracers
308     INTEGER,INTENT(IN) :: nqtot ! # of tracers
309     INTEGER,INTENT(IN) :: ngrid ! # of physics columns
310     REAL,INTENT(OUT) :: pu(ngrid,llm) ! zonal wind speed on physics grid
311     REAL,INTENT(OUT) :: pv(ngrid,llm) ! meridional wind speed on physics grid
312     REAL,INTENT(OUT) :: pr(ngrid,llm) ! vorticity on physics grid
313     REAL,INTENT(OUT) :: pq(ngrid,llm,nqtot) ! tracers on physics grid
314     
315     INTEGER :: i,j,l,iq,ig0
316     REAL :: rot(iim,jjm,llm) ! rot, on dynamics grid
317     REAL :: z1(iim),zsin(iim),zcos(iim)
318     REAL,EXTERNAL :: ssum
319     
320     ! 1. Compute zonal wind on physics grid
321     do l=1,llm
322       do j=2,jjm
323         ig0=1+(j-2)*iim
324         pu(ig0+1,l)=0.5*(ucov(iim,j,l)/cu(iim,j) + ucov(1,j,l)/cu(1,j))
325         do i=2,iim
326           pu(ig0+i,l)=0.5*(ucov(i-1,j,l)/cu(i-1,j) + ucov(i,j,l)/cu(i,j))
327         enddo
328       enddo ! of do j=2,jjm
329     enddo ! of do l=1,llm
330     
331     ! North pole: u = 1/pi * Integral[v * cos(long) * d long]
332     do l=1,llm
333       z1(1)=(rlonu(1)-rlonu(iim)+2.*pi)*vcov(1,1,l)/cv(1,1)
334       do i=2,iim
335         z1(i)=(rlonu(i)-rlonu(i-1))*vcov(i,1,l)/cv(i,1)
336       enddo
337       zcos(1:iim)=cos(rlonv(1:iim))*z1(1:iim)
338       pu(1,l)=ssum(iim,zcos,1)/pi
339     enddo ! of l=1,llm
340     
341     ! South pole: u = 1/pi * Integral[v * cos(long) * d long]
342     do l=1,llm
343       z1(1)=(rlonu(1)-rlonu(iim)+2.*pi)*vcov(1,jjm,l)/cv(1,jjm)
344       do i=2,iim
345         z1(i)=(rlonu(i)-rlonu(i-1))*vcov(i,jjm,l)/cv(i,jjm)
346       enddo
347       zcos(1:iim)=cos(rlonv(1:iim))*z1(1:iim)
348       pu(ngrid,l)=ssum(iim,zcos,1)/pi     
349     enddo ! of do l=1,llm
350     
351     
352     ! 2. Compute meridional wind on physics grid
353     do l=1,llm
354       do j=2,jjm
355         ig0=1+(j-2)*iim
356         do i=1,iim
357           pv(ig0+i,l)=0.5*(vcov(i,j-1,l)/cv(i,j-1) + vcov(i,j,l)/cv(i,j))
358         enddo
359       enddo ! of do j=2,jjm
360     enddo ! of do l=1,llm
361     
362     ! North pole: v = 1/pi * Integral[v * sin(long) * d long ]
363     do l=1,llm
364       z1(1)=(rlonu(1)-rlonu(iim)+2.*pi)*vcov(1,1,l)/cv(1,1)
365       do i=2,iim
366         z1(i)=(rlonu(i)-rlonu(i-1))*vcov(i,1,l)/cv(i,1)
367       enddo
368       zsin(1:iim)=sin(rlonv(1:iim))*z1(1:iim)
369       pv(1,l)=ssum(iim,zsin,1)/pi
370     enddo ! of l=1,llm
371     
372     ! South pole: v = 1/pi * Integral[v * sin(long) * d long ]
373     do l=1,llm
374       z1(1)=(rlonu(1)-rlonu(iim)+2.*pi)*vcov(1,jjm,l)/cv(1,jjm)
375       do i=2,iim
376         z1(i)=(rlonu(i)-rlonu(i-1))*vcov(i,jjm,l)/cv(i,jjm)
377       enddo
378       zsin(1:iim)=sin(rlonv(1:iim))*z1(1:iim)
379       pv(ngrid,l)=ssum(iim,zsin,1)/pi
380     enddo  ! of do l=1,llm
381     
382     
383     ! 3. Compute vorticity, first on dynamics grid
384     do l=1,llm
385       do i=1,iim
386         do j=1,jjm
387           rot(i,j,l)=(vcov(i+1,j,l)-vcov(i,j,l) + ucov(i,j+1,l)-ucov(i,j,l)) / &
388                      ((cu(i,j)+cu(i,j+1))*(cv(i+1,j)+cv(i,j))*4)
389         enddo
390       enddo ! of do i=1,iim
391     enddo  ! of do l=1,llm
392     
393     ! compute vorticity, on physics grid
394     do l=1,llm
395       ! North Pole:
396       pr(1,l)=0.
397       do j=2,jjm
398         ig0=1+(j-2)*iim
399         pr(ig0+1,l)=0.25*(rot(iim,j-1,l)+rot(iim,j,l)+ &
400                           rot(1,j-1,l)+rot(1,j,l))
401         do i=2,iim
402           pr(ig0+i,l)=0.25*(rot(i-1,j-1,l)+rot(i-1,j,l)+ &
403                             rot(i,j-1,l)+rot(i,j,l))
404         enddo
405       enddo ! of do j=2,jjm
406       ! South pole
407       pr(ngrid,l)=0.
408     enddo ! of do l=1,llm
409     
410     ! 4. Tracers, copy them over to the physics grid
411     do iq=1,nqtot
412       CALL gr_dyn_fi(llm,iip1,jjp1,ngrid,q(1,1,1,iq),pq(1,1,iq))
413     enddo
414     
415  END SUBROUTINE init_winds_tracers
416
417!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
418
419  SUBROUTINE init_temperature_pressure_geopot(ps,masse,teta,phis, &
420                                              ngrid,pplev,pplay,pphi,pphis,pt)
421
422    USE comvert_mod, ONLY: ap,bp
423    USE exner_hyb_m, ONLY: exner_hyb
424    USE comconst_mod, ONLY: cpp
425
426    IMPLICIT NONE
427
428    include "dimensions.h"
429    include "paramet.h"
430
431    REAL,INTENT(IN) :: ps(ip1jmp1)   ! surface pressure (Pa) on dyn grid
432    REAL,INTENT(IN) :: masse(ip1jmp1,llm) ! air mass on dyn. grid
433    REAL,INTENT(IN) :: teta(ip1jmp1,llm) ! potential temperature on dyn grid
434    REAL,INTENT(IN) :: phis(ip1jmp1)   ! surface geopotential on dyn grid
435    INTEGER,INTENT(IN) :: ngrid ! # of physics columns
436    REAL,INTENT(OUT) :: pplev(ngrid,llm+1) ! pressure at interfaces
437    REAL,INTENT(OUT) :: pplay(ngrid, llm) ! pressure at mid-layer
438    REAL,INTENT(OUT) :: pphi(ngrid, llm) ! geopotential phi=gz at mid-layer
439    REAL,INTENT(OUT) :: pphis(ngrid) ! surface geopotential
440    REAL,INTENT(OUT) :: pt(ngrid, llm) ! temperature at mid-layer
441
442    REAL :: p(ip1jmp1,llmp1) ! interlayer pressure on dynamics grid
443    REAL :: pks(ip1jmp1) ! Exner at the surface
444    REAL :: pk(ip1jmp1,llm) ! Exner at mid-layer
445    REAL :: phi(ip1jmp1,llm) ! Geopotential on dynamics grid
446   
447    REAL :: ppk(ngrid,llm) ! Exner at mid-layer on physics grid
448    REAL :: pteta(ngrid,llm) ! potential temperature on physics grid
449   
450    INTEGER :: ig,l
451
452    ! 1. Compute necessary intermediate variables on dynamics grid
453    ! Compute pressure on dynamics grid
454    CALL pression(ip1jmp1,ap,bp,ps,p)
455    ! Compute Exner on dynamics grid
456    CALL exner_hyb(ip1jmp1,ps,p,pks,pk)
457    ! Compute geopotential on dynamics grid
458    CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
459   
460    ! 2. Output fields, on physics grid
461   
462    ! Copy over p() to pplev()
463    CALL gr_dyn_fi(llm+1,iip1,jjp1,ngrid,p,pplev)
464   
465    ! Compute pplay() based on ab(),bp() and pplev(:,1) (surface pressure)
466    do ig=1,ngrid
467      pplay(ig,1:llm)=0.5*(ap(1:llm)+bp(1:llm)*pplev(ig,1)+&
468                           ap(2:llm+1)+bp(2:llm+1)*pplev(ig,1))
469    enddo
470   
471    ! Copy surface geopotential phis() to pphis()
472    CALL gr_dyn_fi(1,iip1,jjp1,ngrid,phis,pphis)
473   
474    ! Copy geopotential phi() to pphi()
475    CALL gr_dyn_fi(llm,iip1,jjp1,ngrid,phi,pphi)
476    ! and make pphi relative to surface
477    do l=1,llm
478      pphi(1:ngrid,l)=pphi(1:ngrid,l)-pphis(1:ngrid)
479    enddo
480   
481    ! Compute temperature from teta, on physics grid
482    CALL gr_dyn_fi(llm,iip1,jjp1,ngrid,pk,ppk)
483    CALL gr_dyn_fi(llm,iip1,jjp1,ngrid,teta,pteta)
484    do l=1,llm
485      pt(1:ngrid,l)=pteta(1:ngrid,l)*ppk(1:ngrid,l)/cpp
486    enddo
487   
488
489  END SUBROUTINE init_temperature_pressure_geopot
490
491END PROGRAM driver
Note: See TracBrowser for help on using the repository browser.