source: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90 @ 3807

Last change on this file since 3807 was 3741, checked in by jbclement, 2 months ago

Mars PCM:
Simplification and clarification for the prescription of the atmospheric water profile in 1D: the boolean 'prescribed_h2ovap' activates the option, the file "profile_prescribed_h2ovap" defines the prescribed profile and the value 'h2ovap_relax_time', if positive, activates the relaxation and gives the constant of time.
JBC

File size: 11.3 KB
Line 
1PROGRAM testphys1d
2
3use comsoil_h,           only: inertiedat, inertiesoil, nsoilmx, tsoil, nqsoil, qsoil
4use surfdat_h,           only: albedodat, perennial_co2ice, watercap, tsurf, emis, qsurf
5use comslope_mod,        only: def_slope, subslope_dist
6use phyredem,            only: physdem0, physdem1
7use watersat_mod,        only: watersat
8use tracer_mod,          only: igcm_h2o_vap, igcm_h2o_ice, noms
9use comcstfi_h,          only: pi, g, rcp, cpp
10use time_phylmdz_mod,    only: daysec
11use dimradmars_mod,      only: tauvis, totcloudfrac, albedo
12use dust_param_mod,      only: tauscaling
13use comvert_mod,         only: ap, bp, aps, bps
14use physiq_mod,          only: physiq
15use turb_mod,            only: q2
16use write_output_mod,    only: write_output
17use ioipsl_getincom,     only: getin ! To use 'getin'
18use init_testphys1d_mod, only: init_testphys1d
19use writerestart1D_mod,  only: writerestart1D
20! Mostly for XIOS outputs:
21use mod_const_mpi,       only: init_const_mpi
22use parallel_lmdz,       only: init_parallel
23use version_info_mod,    only: print_version_info
24use callkeys_mod, only: water
25
26implicit none
27
28!=======================================================================
29!   subject:
30!   --------
31!   PROGRAM useful to run physical part of the martian GCM in a 1D column
32!
33! Can be compiled with a command like (e.g. for 25 layers)
34!   "makegcm -p mars -d 25 testphys1d"
35! It requires the files "testphys1d.def" "callphys.def"
36!   and a 'run.def' file (containing a "INCLUDEDEF=callphys.def" line)
37!   and a file describing the sigma layers (e.g. "z2sig.def")
38!
39!   author: Frederic Hourdin, R. Fournier, F. Forget
40!   -------
41!
42!   update: 12/06/2003, including chemistry (S. Lebonnois)
43!                       and water ice (F. Montmessin)
44!           27/09/2023, upgrade to F90 (JB Clément)
45!
46!=======================================================================
47
48include "dimensions.h"
49
50!--------------------------------------------------------------
51! Declarations
52!--------------------------------------------------------------
53integer, parameter                  :: ngrid = 1     ! (2+(jjm-1)*iim - 1/jjm)
54integer, parameter                  :: nlayer = llm
55real, parameter                     :: odpref = 610. ! DOD reference pressure (Pa)
56integer                             :: ndt, ilayer, idt
57logical                             :: firstcall, lastcall
58integer                             :: day0          ! initial (sol ; =0 at Ls=0)
59real                                :: day           ! date during the run
60real                                :: time          ! time (0<time<1 ; time=0.5 a midi)
61real                                :: dttestphys    ! testphys1d timestep
62real, dimension(nlayer)             :: play          ! Pressure at the middle of the layers (Pa)
63real, dimension(nlayer + 1)         :: plev          ! intermediate pressure levels (Pa)
64real                                :: psurf         ! Surface pressure
65real, dimension(nlayer)             :: u, v          ! zonal, meridional wind
66real                                :: gru, grv      ! prescribed "geostrophic" background wind
67real, dimension(nlayer)             :: temp          ! temperature at the middle of the layers
68real, dimension(:,:,:), allocatable :: q             ! tracer mixing ratio (e.g. kg/kg)
69real, dimension(1)                  :: wstar = 0.    ! Thermals vertical velocity
70character(100)                      :: arg           ! To read command-line arguments
71
72! Physical and dynamical tendencies (e.g. m.s-2, K/s, Pa/s)
73real, dimension(nlayer)             :: du, dv, dtemp
74real, dimension(1)                  :: dpsurf
75real, dimension(:,:,:), allocatable :: dq, dqdyn
76
77! Various intermediate variables
78real                    :: zls, pks, ptif
79real, dimension(nlayer) :: phi, h, s, w
80integer                 :: nq = 1 ! number of tracers
81real, dimension(1)      :: latitude, longitude, cell_area
82
83! Use of starting files for 1D
84logical :: startfiles_1D, therestart1D, therestartfi, there
85
86! Prescription of atmospheric water profile
87logical                         :: prescribed_h2ovap
88real                            :: h2ovap_relax_time
89real, dimension(nlayer)         :: q_prescribed_h2o_vap
90real, dimension(:), allocatable :: zqsat
91!=======================================================================
92
93!=======================================================================
94! INITIALISATION
95!=======================================================================
96if (command_argument_count() > 0) then ! Get the number of command-line arguments
97    call get_command_argument(1,arg) ! Read the argument given to the program
98    select case (trim(adjustl(arg)))
99        case('version')
100            call print_version_info()
101            stop
102        case default
103            error stop 'The argument given to the program is unknown!'
104    end select
105endif
106
107#ifdef CPP_XIOS
108    call init_const_mpi
109    call init_parallel
110#endif
111
112! Initialize "serial/parallel" related stuff
113!call init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
114!call init_phys_lmdz(1,1,llm,1,(/1/))
115!call initcomgeomphy
116
117startfiles_1D = .false.
118!------------------------------------------------------
119! Loading run parameters from "run.def" file
120!------------------------------------------------------
121! check if 'run.def' file is around. Otherwise reading parameters
122! from callphys.def via getin() routine won't work.
123inquire(file = 'run.def',exist = there)
124if (.not. there) then
125    write(*,*) 'Cannot find required file "run.def"'
126    write(*,*) '  (which should contain some input parameters along with the following line: INCLUDEDEF=callphys.def)'
127    write(*,*) ' ... might as well stop here ...'
128    error stop
129endif
130
131startfiles_1D = .false. ! Default value
132write(*,*) 'Do you want to use starting files and/or to write restarting files?'
133call getin('startfiles_1D',startfiles_1D)
134write(*,*) 'startfiles_1D =', startfiles_1D
135
136therestart1D = .false. ! Default value
137therestartfi = .false. ! Default value
138if (startfiles_1D) then
139    inquire(file = 'start1D.txt',exist = therestart1D)
140    if (startfiles_1D .and. .not. therestart1D) then
141        write(*,*) 'There is no "start1D.txt" file!'
142        write(*,*) 'Initialization is done with default values.'
143    endif
144    inquire(file = 'startfi.nc',exist = therestartfi)
145    if (.not. therestartfi) then
146        write(*,*) 'There is no "startfi.nc" file!'
147        write(*,*) 'Initialization is done with default values.'
148    endif
149endif
150
151call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,odpref,nq,q, &
152                     time,psurf,u,v,temp,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
153                     play,plev,latitude,longitude,cell_area,prescribed_h2ovap,h2ovap_relax_time,q_prescribed_h2o_vap)
154
155! Write a "startfi" file
156! ----------------------
157! This file will be read during the first call to "physiq".
158! It is needed to transfert physics variables to "physiq"...
159if (.not. therestartfi) then
160    call physdem0("startfi.nc",longitude,latitude,nsoilmx,ngrid, &
161                  llm,nq,dttestphys,float(day0),0.,cell_area,    &
162                  albedodat,inertiedat,def_slope,subslope_dist)
163    call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,nqsoil,dttestphys,time,      &
164                  tsurf,tsoil,inertiesoil,albedo,emis,q2,qsurf,qsoil,tauscaling, &
165                  totcloudfrac,wstar,watercap,perennial_co2ice)
166endif !(.not. therestartfi)
167
168!=======================================================================
169!  1D MODEL TIME STEPPING LOOP
170!=======================================================================
171firstcall = .true.
172lastcall = .false.
173do idt = 1,ndt
174    if (idt == ndt) lastcall = .true.
175!    if (idt == ndt - day_step - 1) then ! test
176!        lastcall = .true.
177!        call solarlong(day*1.0,zls)
178!        write(103,*) 'Ls=',zls*180./pi
179!        write(103,*) 'Lat=', latitude(1)*180./pi
180!        write(103,*) 'Tau=', tauvis/odpref*psurf
181!        write(103,*) 'RunEnd - Atmos. Temp. File'
182!        write(103,*) 'RunEnd - Atmos. Temp. File'
183!        write(104,*) 'Ls=',zls*180./pi
184!        write(104,*) 'Lat=', latitude(1)
185!        write(104,*) 'Tau=', tauvis/odpref*psurf
186!        write(104,*) 'RunEnd - Atmos. Temp. File'
187!    endif
188
189    ! Compute geopotential
190    ! ~~~~~~~~~~~~~~~~~~~~
191    s = (aps/psurf + bps)**rcp
192    h = cpp*temp/(pks*s)
193
194    phi(1) = pks*h(1)*(1. - s(1))
195    do ilayer = 2,nlayer
196        phi(ilayer) = phi(ilayer - 1) + pks*(h(ilayer - 1) + h(ilayer))*.5*(s(ilayer - 1) - s(ilayer))
197    enddo
198
199    ! Prescription of atmospheric water if asked
200    ! ------------------------------------------
201    if (water) then
202        call watersat(nlayer,temp,play,zqsat)
203        if (prescribed_h2ovap) then ! If atmospheric water profile is prescribed
204            if (h2ovap_relax_time < 0.) then ! Forcing
205                q(1,:,igcm_h2o_vap) = min(zqsat(:),q_prescribed_h2o_vap(:))
206            else ! Relaxation
207                q(1,:,igcm_h2o_vap) = q_prescribed_h2o_vap(:) + (q(1,:,igcm_h2o_vap) - q_prescribed_h2o_vap(:))*dexp(-dttestphys/h2ovap_relax_time)
208            endif
209        endif
210     endif
211
212    ! Call physics
213    ! --------------------
214    call physiq(1,llm,nq,firstcall,lastcall,day,time,dttestphys,plev,play,phi,u,v,temp,q,w,du,dv,dtemp,dq,dpsurf)
215    !                                                                                      ^----- outputs -----^
216
217    ! Wind increment: specific for 1D
218    ! --------------------------------
219    ! The physics compute the tendencies on u and v,
220    ! here we just add Coriolos effect
221    !do ilayer = 1,nlayer
222    !    du(ilayer) = du(ilayer) + ptif*(v(ilayer) - grv)
223    !    dv(ilayer) = dv(ilayer) + ptif*(-u(ilayer) + gru)
224    !enddo
225    ! For some tests: No coriolis force at equator
226    !if (latitude(1) == 0.) then
227    du = du + (gru - u)/1.e4
228    dv = dv + (grv - v)/1.e4
229    !endif
230
231    ! Compute time for next time step
232    ! -------------------------------
233    firstcall = .false.
234    time = time + dttestphys/daysec
235    if (time > 1.) then
236        time = time - 1.
237        day = day + 1
238    endif
239
240    ! Compute winds and temperature for next time step
241    ! ------------------------------------------------
242    u = u + dttestphys*du
243    v = v + dttestphys*dv
244    temp = temp + dttestphys*dtemp
245
246    ! Compute pressure for next time step
247    ! -----------------------------------
248    psurf = psurf + dttestphys*dpsurf(1) ! Surface pressure change
249    plev = ap + psurf*bp
250    play = aps + psurf*bps
251
252    ! Increment tracers
253    q = q + dttestphys*dq
254enddo ! End of time stepping loop (idt=1,ndt)
255
256! Writing the "restart1D.txt" file for the next run
257if (startfiles_1D) call writerestart1D('restart1D.txt',psurf,tsurf(1,:),nlayer,size(tsurf,2),temp,u,v,nq,noms,qsurf(1,:,:),q)
258
259write(*,*) "testphys1d: everything is cool!"
260
261END PROGRAM testphys1d
262
263!***********************************************************************
264!***********************************************************************
265! Dummy subroutine used only in 3D, but required to
266! compile testphys1d (to cleanly use writediagfi)
267SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
268
269implicit none
270
271integer                       :: im, jm, ngrid, nfield
272real, dimension(im,jm,nfield) :: pdyn
273real, dimension(ngrid,nfield) :: pfi
274
275if (ngrid /= 1) error stop 'gr_fi_dyn error: in 1D ngrid should be 1!!!'
276
277pdyn(1,1,1:nfield) = pfi(1,1:nfield)
278
279END SUBROUTINE gr_fi_dyn
280
281!***********************************************************************
282!***********************************************************************
Note: See TracBrowser for help on using the repository browser.