source: trunk/LMDZ.MARS/libf/phymars/phyredem.F90 @ 2266

Last change on this file since 2266 was 2265, checked in by emillour, 5 years ago

Mars GCM:
Save "dtau", the opacity difference between model and target dust scenario
in the restartfi.nc file so that we have 1+1=2 when running with dust
injection schemes.
EM

File size: 10.3 KB
Line 
1module phyredem
2
3implicit none
4
5contains
6
7subroutine physdem0(filename,lonfi,latfi,nsoil,ngrid,nlay,nq, &
8                         phystep,day_ini,time,airefi, &
9                         alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe, &
10                         phmons,psummit,pbase)
11! create physics restart file and write time-independent variables
12  use comsoil_h, only: inertiedat, volcapa, mlayer
13  use geometry_mod, only: cell_area
14  use surfdat_h, only: zmea, zstd, zsig, zgam, zthe, &
15                       z0_default, albedice, emisice, emissiv, &
16                       iceradius, dtemisice, phisfi, z0,   &
17                       hmons,summit,base
18  use dimradmars_mod, only: tauvis
19  use iostart, only : open_restartphy, close_restartphy, &
20                      put_var, put_field, length
21  use mod_grid_phy_lmdz, only : klon_glo
22  use planete_h, only: aphelie, emin_turb, lmixmin, obliquit, &
23                       peri_day, periheli, year_day
24  use comcstfi_h, only: g, mugaz, omeg, rad, rcp
25  use time_phylmdz_mod, only: daysec
26  implicit none
27 
28  character(len=*), intent(in) :: filename
29  real,intent(in) :: lonfi(ngrid)
30  real,intent(in) :: latfi(ngrid)
31  integer,intent(in) :: nsoil
32  integer,intent(in) :: ngrid
33  integer,intent(in) :: nlay
34  integer,intent(in) :: nq
35  real,intent(in) :: phystep
36  real,intent(in) :: day_ini
37  real,intent(in) :: time
38  real,intent(in) :: airefi(ngrid)
39  real,intent(in) :: alb(ngrid)
40  real,intent(in) :: ith(ngrid,nsoil)
41  real,intent(in) :: pzmea(ngrid)
42  real,intent(in) :: pzstd(ngrid)
43  real,intent(in) :: pzsig(ngrid)
44  real,intent(in) :: pzgam(ngrid)
45  real,intent(in) :: pzthe(ngrid)
46  real,intent(in) :: phmons(ngrid)
47  real,intent(in) :: psummit(ngrid)
48  real,intent(in) :: pbase(ngrid)
49
50  real :: tab_cntrl(length) ! nb "length=100" defined in iostart module
51 
52  ! Create physics start file
53  call open_restartphy(filename)
54 
55  ! Build tab_cntrl(:) array
56  tab_cntrl(:)=0.0
57  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
58  ! Fill control array tab_cntrl(:) with parameters for this run
59  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
60  ! Informations on the physics grid
61  tab_cntrl(1) = float(klon_glo)  ! total number of nodes on physics grid
62  tab_cntrl(2) = float(nlay) ! number of atmospheric layers
63  tab_cntrl(3) = day_ini + int(time)      ! initial day
64  tab_cntrl(4) = time -int(time)          ! initial time of day
65
66  ! Informations about Mars, used by dynamics and physics
67  tab_cntrl(5) = rad      ! radius of Mars (m) ~3397200
68  tab_cntrl(6) = omeg     ! rotation rate (rad.s-1)
69  tab_cntrl(7) = g        ! gravity (m.s-2) ~3.72
70  tab_cntrl(8) = mugaz    ! Molar mass of the atmosphere (g.mol-1) ~43.49
71  tab_cntrl(9) = rcp      !  = r/cp  ~0.256793 (=kappa dans dynamique)
72  tab_cntrl(10) = daysec  ! length of a sol (s)  ~88775
73
74  tab_cntrl(11) = phystep  ! time step in the physics
75  tab_cntrl(12) = 0.
76  tab_cntrl(13) = 0.
77
78  ! Informations about Mars, only for physics
79  tab_cntrl(14) = year_day  ! length of year (sols) ~668.6
80  tab_cntrl(15) = periheli  ! min. Sun-Mars distance (Mkm) ~206.66
81  tab_cntrl(16) = aphelie   ! max. SUn-Mars distance (Mkm) ~249.22
82  tab_cntrl(17) = peri_day  ! date of perihelion (sols since N. spring)
83  tab_cntrl(18) = obliquit  ! Obliquity of the planet (deg) ~23.98
84
85  ! Boundary layer and turbulence
86  tab_cntrl(19) = z0_default   ! default surface roughness (m) ~0.01
87  tab_cntrl(20) = lmixmin   ! mixing length ~100
88  tab_cntrl(21) = emin_turb ! minimal energy ~1.e-8
89
90  ! Optical properties of polar caps and ground emissivity
91  tab_cntrl(22) = albedice(1)  ! Albedo of northern cap ~0.5
92  tab_cntrl(23) = albedice(2)  ! Albedo of southern cap ~0.5
93  tab_cntrl(24) = emisice(1)   ! Emissivity of northern cap ~0.95
94  tab_cntrl(25) = emisice(2)   ! Emissivity of southern cap ~0.95
95  tab_cntrl(26) = emissiv      ! Emissivity of martian soil ~.95
96  tab_cntrl(31) = iceradius(1) ! mean scat radius of CO2 snow (north)
97  tab_cntrl(32) = iceradius(2) ! mean scat radius of CO2 snow (south)
98  tab_cntrl(33) = dtemisice(1) ! time scale for snow metamorphism (north)
99  tab_cntrl(34) = dtemisice(2) ! time scale for snow metamorphism (south)
100
101  ! dust aerosol properties
102  tab_cntrl(27) = tauvis      ! mean visible optical depth
103
104  ! Soil properties:
105  tab_cntrl(35) = volcapa ! soil volumetric heat capacity
106
107  ! Write the controle array
108  call put_var("controle","Control parameters",tab_cntrl)
109 
110  ! Write the mid-layer depths
111  call put_var("soildepth","Soil mid-layer depth",mlayer)
112 
113  ! Write longitudes
114  call put_field("longitude","Longitudes of physics grid",lonfi)
115 
116  ! Write latitudes
117  call put_field("latitude","Latitudes of physics grid",latfi)
118 
119  ! Write mesh areas
120  call put_field("area","Mesh area",cell_area)
121 
122  ! Write surface geopotential
123  call put_field("phisfi","Geopotential at the surface",phisfi)
124 
125  ! Write surface albedo
126  call put_field("albedodat","Albedo of bare ground",alb)
127 
128  ! Subgrid topogaphy variables
129  call put_field("ZMEA","Relief: mean relief",zmea)
130  call put_field("ZSTD","Relief: standard deviation",zstd)
131  call put_field("ZSIG","Relief: sigma parameter",zsig)
132  call put_field("ZGAM","Relief: gamma parameter",zgam)
133  call put_field("ZTHE","Relief: theta parameter",zthe)
134  call put_field("hmons","Relief: hmons parameter (summit - base)",hmons)
135  call put_field("summit","Relief: altitude from the aeroid of the summit of the highest subgrid topography",summit)
136  call put_field("base","Relief: altitude from the aeroid of the base of the highest subgrid topography",base)
137     
138  ! Soil thermal inertia
139  call put_field("inertiedat","Soil thermal inertia",ith)
140 
141  ! Surface roughness length
142  call put_field("z0","Surface roughness length",z0)
143 
144  ! Close file
145  call close_restartphy
146 
147end subroutine physdem0
148
149subroutine physdem1(filename,nsoil,ngrid,nlay,nq, &
150                    phystep,time,tsurf,tsoil,co2ice,albedo,emis,q2,qsurf,&
151                    tauscaling,totcloudfrac,wstar, &
152                    mem_Mccn_co2,mem_Nccn_co2,mem_Mh2o_co2, watercap)
153  ! write time-dependent variable to restart file
154  use iostart, only : open_restartphy, close_restartphy, &
155                      put_var, put_field
156  use tracer_mod, only: noms ! tracer names
157  use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd
158  use compute_dtau_mod, only: dtau
159
160  implicit none
161 
162  include "callkeys.h"
163 
164  character(len=*),intent(in) :: filename
165  integer,intent(in) :: nsoil
166  integer,intent(in) :: ngrid
167  integer,intent(in) :: nlay
168  integer,intent(in) :: nq
169  real,intent(in) :: phystep
170  real,intent(in) :: time
171  real,intent(in) :: tsurf(ngrid)
172  real,intent(in) :: tsoil(ngrid,nsoil)
173  real,intent(in) :: co2ice(ngrid)
174  real,intent(in) :: albedo(ngrid,2)
175  real,intent(in) :: emis(ngrid)
176  real,intent(in) :: q2(ngrid,nlay+1)
177  real,intent(in) :: qsurf(ngrid,nq)
178  real,intent(in) :: tauscaling(ngrid)
179  real,intent(in) :: totcloudfrac(ngrid)
180  real,intent(in) :: wstar(ngrid)
181  real,intent(in) :: mem_Mccn_co2(ngrid,nlay) ! CCN mass of H2O and dust used by CO2
182  real,intent(in) :: mem_Nccn_co2(ngrid,nlay) ! CCN number of H2O and dust used by CO2
183  real,intent(in) :: mem_Mh2o_co2(ngrid,nlay) ! H2O mass integred into CO2 crystal
184  real,intent(in) :: watercap(ngrid)
185 
186  integer :: iq
187  character(len=30) :: txt ! to store some text
188  ! indexes of water vapour & water ice tracers (if any):
189  integer :: i_h2o_vap=0
190  integer :: i_h2o_ice=0
191
192 
193  ! Open file
194  call open_restartphy(filename)
195 
196  ! First variable to write must be "Time", in order to correctly
197  ! set time counter in file
198  call put_var("Time","Temps de simulation",time)
199 
200  ! CO2 ice layer
201  call put_field("co2ice","CO2 ice cover",co2ice,time)
202
203  ! Water ice layer
204  call put_field("watercap","H2O ice cover",watercap,time)
205 
206  ! Surface temperature
207  call put_field("tsurf","Surface temperature",tsurf,time)
208 
209  ! Soil temperature
210  call put_field("tsoil","Soil temperature",tsoil,time)
211 
212  ! Albedo of the surface
213  call put_field("albedo","Surface albedo",albedo(:,1),time)
214 
215  ! Emissivity of the surface
216  call put_field("emis","Surface emissivity",emis,time)
217 
218  ! Planetary Boundary Layer
219  call put_field("q2","pbl wind variance",q2,time)
220
221  ! Sub-grid cloud fraction
222  call put_field("totcloudfrac","Total cloud fraction",totcloudfrac,time)
223 
224  ! Dust conversion factor
225  ! Only to be read by newstart to convert to actual dust values
226  ! Or by any user who wants to reconstruct dust, opacity from the start files.
227  call put_field("tauscaling","dust conversion factor",tauscaling,time)
228
229  if (dustinjection.gt.0) then
230    call put_field("dtau","dust opacity difference between GCM and scenario",&
231                   dtau,time)
232  endif
233
234  if (calltherm) then
235    call put_field("wstar","Max vertical velocity in thermals",wstar,time)
236  endif
237
238  ! Tracers on the surface
239  ! preliminary stuff: look for water vapour & water ice tracers (if any)
240  do iq=1,nq
241    if (noms(iq).eq."h2o_vap") then
242      i_h2o_vap=iq
243    endif
244    if (noms(iq).eq."h2o_ice") then
245      i_h2o_ice=iq
246    endif
247  enddo
248 
249  if (nq.gt.0) then
250    do iq=1,nq
251      txt=noms(iq)
252      ! Exception: there is no water vapour surface tracer
253      if (txt.eq."h2o_vap") then
254        write(*,*)"physdem1: skipping water vapour tracer"
255        if (i_h2o_ice.eq.i_h2o_vap) then
256          ! then there is no "water ice" tracer; but still
257          ! there is some water ice on the surface
258          write(*,*)"          writing water ice instead"
259          txt="h2o_ice"
260        else
261          ! there is a "water ice" tracer which has been / will be
262          ! delt with in due time
263          cycle
264        endif ! of if (igcm_h2o_ice.eq.igcm_h2o_vap)
265      endif ! of if (txt.eq."h2o_vap")
266      call put_field(trim(txt),"tracer on surface",qsurf(:,iq),time)
267    enddo
268  endif
269  ! Memory of the origin of the co2 particles
270  if (co2useh2o) then
271     call put_field("mem_Mccn_co2","CCN mass of H2O and dust used by CO2",mem_Mccn_co2,time)
272     call put_field("mem_Nccn_co2","CCN number of H2O and dust used by CO2",mem_Nccn_co2,time)
273     call put_field("mem_Mh2o_co2","H2O mass integred into CO2 crystal",mem_Mh2o_co2,time)
274  endif
275 
276  ! Non-orographic gavity waves
277  if (calllott_nonoro) then
278     call put_field("du_nonoro_gwd","Zonal wind tendency due to GW",du_nonoro_gwd,time)
279     call put_field("dv_nonoro_gwd","Meridional wind tendency due to GW",dv_nonoro_gwd,time)
280  endif
281  ! Close file
282  call close_restartphy
283 
284end subroutine physdem1
285
286end module phyredem
Note: See TracBrowser for help on using the repository browser.