source: trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90 @ 3094

Last change on this file since 3094 was 3094, checked in by emillour, 14 months ago

Mars PCM:
Fix issues with dates when in parallel with OpenMP (missing copyin
statement when opening parallel section in iniphysiq). Added some threadprivate
clauses in saved module variables in comcstfi_h.F90 and planete_h.F90.
Prettyfied solarlong.F and made it a module. Likewise for conf_phys.F
EM

File size: 28.5 KB
Line 
1MODULE init_testphys1d_mod
2
3implicit none
4
5contains
6
7SUBROUTINE init_testphys1d(pem1d,ngrid,nlayer,odpref,nq,q,time,psurf,u,v,temp,startfiles_1D,therestart1D, &
8                           therestartfi,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,        &
9                           play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
10
11use ioipsl_getincom,          only: getin ! To use 'getin'
12use comcstfi_h,               only: pi, rad, omeg, g, mugaz, rcp, r, cpp
13use time_phylmdz_mod,         only: daysec, day_step, ecritphy, iphysiq
14use planete_h,                only: year_day, periheli, aphelie, peri_day, obliquit, emin_turb, lmixmin
15use surfdat_h,                only: albedodat, z0_default, z0, emissiv, emisice, albedice, iceradius, dtemisice,      &
16                                    zmea, zstd, zsig, zgam, zthe, hmons, summit, base, phisfi, watercaptag, watercap, &
17                                    tsurf, emis, qsurf
18use infotrac,                 only: nqtot, tname, nqperes, nqfils
19use read_profile_mod,         only: read_profile
20use iostart,                  only: open_startphy, get_var, close_startphy
21use physics_distribution_mod, only: init_physics_distribution
22use comsoil_h,                only: volcapa, nsoilmx, inertiesoil, inertiedat, layer, mlayer, flux_geo, tsoil
23use comvert_mod,              only: ap, bp, aps, bps, pa, preff, presnivs, pseudoalt, scaleheight
24use dimradmars_mod,           only: tauvis, totcloudfrac, albedo
25use regular_lonlat_mod,       only: init_regular_lonlat
26use mod_interface_dyn_phys,   only: init_interface_dyn_phys
27use geometry_mod,             only: init_geometry, init_geometry_cell_area_for_outputs
28use dimphy,                   only: init_dimphy
29use comgeomfi_h,              only: sinlat, ini_fillgeom
30use slope_mod,                only: theta_sl, psi_sl
31use comslope_mod,             only: def_slope, subslope_dist
32use dust_param_mod,           only: tauscaling
33use tracer_mod,               only: igcm_co2
34use logic_mod,                only: hybrid
35use vertical_layers_mod,      only: init_vertical_layers
36use inichim_newstart_mod,     only: inichim_newstart
37use mod_grid_phy_lmdz,        only: regular_lonlat
38use phys_state_var_init_mod,  only: phys_state_var_init
39use turb_mod,                 only: q2
40use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd
41use conf_phys_mod, only: conf_phys
42! Mostly for XIOS outputs:
43use mod_const_mpi,            only: COMM_LMDZ
44
45implicit none
46
47include "dimensions.h"
48include "callkeys.h"
49
50!=======================================================================
51! Arguments
52!=======================================================================
53integer, intent(in) :: ngrid, nlayer
54real,    intent(in) :: odpref        ! DOD reference pressure (Pa)
55logical, intent(in) :: pem1d         ! If initialization for the 1D PEM
56
57integer, intent(inout) :: nq
58
59real, dimension(:,:,:), allocatable, intent(out) :: q     ! tracer mixing ratio (e.g. kg/kg)
60real,                                intent(out) :: time  ! time (0<time<1; time=0.5 at noon)
61real,                                intent(out) :: psurf ! Surface pressure
62real, dimension(nlayer),             intent(out) :: u, v  ! zonal, meridional wind
63real, dimension(nlayer),             intent(out) :: temp  ! temperature at the middle of the layers
64logical,                             intent(out) :: startfiles_1D, therestart1D, therestartfi ! Use of starting files for 1D
65integer,                             intent(out) :: ndt
66real,                                intent(out) :: ptif, pks
67real,                                intent(out) :: dttestphys                   ! testphys1d timestep
68real, dimension(:),     allocatable, intent(out) :: zqsat                        ! useful for (atm_wat_profile=2)
69real, dimension(:,:,:), allocatable, intent(out) :: dq, dqdyn                    ! Physical and dynamical tandencies
70integer,                             intent(out) :: day0                         ! initial (sol ; =0 at Ls=0) and final date
71real,                                intent(out) :: day                          ! date during the run
72real,                                intent(out) :: gru, grv                     ! prescribed "geostrophic" background wind
73real, dimension(nlayer),             intent(out) :: w                            ! "Dummy wind" in 1D
74real, dimension(nlayer),             intent(out) :: play                         ! Pressure at the middle of the layers (Pa)
75real, dimension(nlayer + 1),         intent(out) :: plev                         ! intermediate pressure levels (pa)
76real, dimension(1),                  intent(out) :: latitude, longitude, cell_area
77real,                                intent(out) :: atm_wat_profile, atm_wat_tau ! Force atmospheric water profiles
78
79!=======================================================================
80! Local variables
81!=======================================================================
82integer                   :: ierr, iq, j, ilayer, isoil, nlevel, nsoil, flagthermo, flagh2o
83integer                   :: dayn ! Final date
84real, dimension(nlayer)   :: zlay ! altitude estimee dans les couches (km)
85real, dimension(0:nlayer) :: tmp1, tmp2
86
87! Dummy variables along "dynamics scalar grid"
88real, dimension(:,:,:,:), allocatable :: qdyn
89real, dimension(:,:),     allocatable :: psdyn
90
91! RV & JBC: Use of starting files for 1D
92logical                :: found
93character(len = 30)    :: header
94real, dimension(100)   :: tab_cntrl
95real, dimension(1,2,1) :: albedo_read ! surface albedo
96
97! New flag to compute paleo orbital configurations + few variables JN
98logical :: paleomars
99real    :: halfaxe, eccentric, Lsperi
100
101! MVals: isotopes as in the dynamics (CRisi)
102integer                                        :: ifils, ipere, generation, ierr0
103character(len = 30), dimension(:), allocatable :: tnom_transp ! transporting fluid short name
104character(len = 80)                            :: line        ! to store a line of text
105logical                                        :: continu, there
106
107! LL: Possibility to add subsurface ice
108real    :: ice_depth          ! depth of the ice table, ice_depth < 0. means no ice
109real    :: inertieice = 2100. ! ice thermal inertia
110integer :: iref
111
112! LL: Subsurface geothermal flux
113real :: flux_geo_tmp
114
115! JBC: To initialize the 1D PEM
116character(:), allocatable :: start1Dname, startfiname ! Name of starting files for 1D
117
118!=======================================================================
119! Code
120!=======================================================================
121if (.not. pem1d) then
122    start1Dname = 'start1D.txt'
123    startfiname = 'startfi.nc'
124    startfiles_1D = .false.
125    !------------------------------------------------------
126    ! Loading run parameters from "run.def" file
127    !------------------------------------------------------
128    ! check if 'run.def' file is around. Otherwise reading parameters
129    ! from callphys.def via getin() routine won't work.
130    inquire(file = 'run.def',exist = there)
131    if (.not. there) then
132        write(*,*) 'Cannot find required file "run.def"'
133        write(*,*) '  (which should contain some input parameters along with the following line: INCLUDEDEF=callphys.def)'
134        write(*,*) ' ... might as well stop here ...'
135        error stop
136    endif
137
138    write(*,*)'Do you want to use starting files?'
139    call getin("startfiles_1D",startfiles_1D)
140    write(*,*) " startfiles_1D = ", startfiles_1D
141else
142    start1dname = 'start1D_evol.txt'
143    startfiname = 'startfi_evol.nc'
144    startfiles_1D = .true.
145endif
146
147therestart1D = .false.
148therestartfi = .false.
149inquire(file = start1Dname,exist = therestart1D)
150if (startfiles_1D .and. .not. therestart1D) then
151    write(*,*) 'There is no "'//start1Dname//'" file!'
152    if (.not. pem1d) then
153        write(*,*) 'Initialization is done with default values.'
154    else
155        error stop 'Initialization cannot be done for the 1D PEM.'
156    endif
157endif
158inquire(file = startfiname,exist = therestartfi)
159if (.not. therestartfi) then
160    write(*,*) 'There is no "'//startfiname//'" file!'
161    if (.not. pem1d) then
162        write(*,*) 'Initialization is done with default values.'
163    else
164        error stop 'Initialization cannot be done for the 1D PEM.'
165    endif
166endif
167
168!------------------------------------------------------
169! Prescribed constants to be set here
170!------------------------------------------------------
171pi = 2.*asin(1.)
172
173! Mars planetary constants
174! ------------------------
175rad = 3397200.            ! mars radius (m) ~3397200 m
176daysec = 88775.           ! length of a sol (s) ~88775 s
177omeg = 4.*asin(1.)/daysec ! rotation rate (rad.s-1)
178g = 3.72                  ! gravity (m.s-2) ~3.72
179mugaz = 43.49             ! atmosphere mola mass (g.mol-1) ~43.49
180rcp = .256793             ! = r/cp ~0.256793
181r = 8.314511*1000./mugaz
182cpp = r/rcp
183year_day = 669            ! length of year (sols) ~668.6
184periheli = 206.66         ! minimum sun-mars distance (Mkm) ~206.66
185aphelie = 249.22          ! maximum sun-mars distance (Mkm) ~249.22
186halfaxe = 227.94          ! demi-grand axe de l'ellipse
187peri_day = 485.           ! perihelion date (sols since N. Spring)
188obliquit = 25.2           ! Obliquity (deg) ~25.2
189eccentric = 0.0934        ! Eccentricity (0.0934)
190
191! Planetary Boundary Layer and Turbulence parameters
192! --------------------------------------------------
193z0_default = 1.e-2 ! surface roughness (m) ~0.01
194emin_turb = 1.e-6  ! minimal turbulent energy ~1.e-8
195lmixmin = 30       ! mixing length ~100
196
197! cap properties and surface emissivities
198! ---------------------------------------
199emissiv = 0.95         ! Bare ground emissivity ~.95
200emisice(1) = 0.95      ! Northern cap emissivity
201emisice(2) = 0.95      ! Southern cap emisssivity
202albedice(1) = 0.5      ! Northern cap albedo
203albedice(2) = 0.5      ! Southern cap albedo
204iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north)
205iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south)
206dtemisice(1) = 2.      ! time scale for snow metamorphism (north)
207dtemisice(2) = 2.      ! time scale for snow metamorphism (south)
208
209! mesh surface (not a very usefull quantity in 1D)
210! ------------------------------------------------
211cell_area(1) = 1.
212
213! check if there is a 'traceur.def' file and process it
214! load tracer names from file 'traceur.def'
215open(90,file = 'traceur.def',status = 'old',form = 'formatted',iostat = ierr)
216if (ierr /= 0) then
217    write(*,*) 'Cannot find required file "traceur.def"'
218    write(*,*) ' If you want to run with tracers, I need it'
219    write(*,*) ' ... might as well stop here ...'
220    error stop
221else
222    write(*,*) "init_testphys1d: Reading file traceur.def"
223    ! read number of tracers:
224    read(90,*,iostat = ierr) nq
225    nqtot = nq ! set value of nqtot (in infotrac module) as nq
226    if (ierr /= 0) then
227        write(*,*) "init_testphys1d: error reading number of tracers"
228        write(*,*) "   (first line of traceur.def) "
229        error stop
230    endif
231    if (nq < 1) then
232        write(*,*) "init_testphys1d: error number of tracers"
233        write(*,*) "is nq=",nq," but must be >=1!"
234        error stop
235    endif
236endif
237
238! allocate arrays:
239allocate(tname(nq),q(1,nlayer,nq),zqsat(nlayer))
240allocate(dq(1,nlayer,nq),dqdyn(1,nlayer,nq),tnom_transp(nq))
241
242! read tracer names from file traceur.def
243do iq = 1,nq
244    read(90,'(80a)',iostat = ierr) line ! store the line from traceur.def
245    if (ierr /= 0) then
246        error stop 'init_testphys1d: error reading tracer names...'
247    endif
248    ! if format is tnom_0, tnom_transp (isotopes)
249    read(line,*,iostat = ierr0) tname(iq),tnom_transp(iq)
250    if (ierr0 /= 0) then
251        read(line,*) tname(iq)
252        tnom_transp(iq) = 'air'
253    endif
254enddo
255close(90)
256
257! Isotopes: as in the 3D case we have to determine father/son relations for isotopes and carrying fluid
258allocate(nqfils(nqtot))
259nqperes = 0
260nqfils(:) = 0
261do iq = 1,nqtot
262    if (tnom_transp(iq) == 'air') then
263        ! ceci est un traceur père
264        write(*,*) 'Le traceur',iq,', appele ',trim(tname(iq)),', est un pere'
265        nqperes = nqperes + 1
266    else !if (tnom_transp(iq) == 'air') then
267        ! ceci est un fils. Qui est son père?
268        write(*,*) 'Le traceur',iq,', appele ',trim(tname(iq)),', est un fils'
269        continu = .true.
270        ipere = 1
271        do while (continu)
272            if (tnom_transp(iq) == tname(ipere)) then
273                ! Son père est ipere
274                write(*,*) 'Le traceur',iq,'appele ',trim(tname(iq)),' est le fils de ',ipere,'appele ',trim(tname(ipere))
275                nqfils(ipere) = nqfils(ipere) + 1
276                continu = .false.
277            else !if (tnom_transp(iq) == tnom_0(ipere)) then
278                ipere = ipere + 1
279                if (ipere > nqtot) then
280                    write(*,*) 'Le traceur',iq,'appele ',trim(tname(iq)),', est orpelin.'
281                    call abort_gcm('infotrac_init','Un traceur est orphelin',1)
282                endif !if (ipere > nqtot) then
283            endif !if (tnom_transp(iq) == tnom_0(ipere)) then
284        enddo !do while (continu)
285    endif !if (tnom_transp(iq) == 'air') then
286enddo !do iq=1,nqtot
287write(*,*) 'nqperes=',nqperes
288write(*,*) 'nqfils=',nqfils
289
290#ifdef CPP_XIOS
291    call init_physics_distribution(regular_lonlat,4,1,1,1,nlayer,COMM_LMDZ)
292#else
293    call init_physics_distribution(regular_lonlat,4,1,1,1,nlayer,1)
294#endif
295
296! Date and local time at beginning of run
297! ---------------------------------------
298if (.not. startfiles_1D) then
299    ! Date (in sols since spring solstice) at beginning of run
300    day0 = 0 ! default value for day0
301    write(*,*) 'Initial date (in martian sols; =0 at Ls=0)?'
302    call getin("day0",day0)
303    day = float(day0)
304    write(*,*) " day0 = ",day0
305    ! Local time at beginning of run
306    time = 0 ! default value for time
307    write(*,*)'Initial local time (in hours, between 0 and 24)?'
308    call getin("time",time)
309    write(*,*)" time = ",time
310    time = time/24. ! convert time (hours) to fraction of sol
311else
312    call open_startphy(startfiname)
313    call get_var("controle",tab_cntrl,found)
314    if (.not. found) then
315        call abort_physic("open_startphy","tabfi: Failed reading <controle> array",1)
316    else
317        write(*,*)'tabfi: tab_cntrl',tab_cntrl
318    endif
319    day0 = tab_cntrl(3)
320    day = float(day0)
321
322    call get_var("Time",time,found)
323    call close_startphy
324endif !startfiles_1D
325
326! Discretization (Definition of grid and time steps)
327! --------------------------------------------------
328nlevel = nlayer + 1
329nsoil = nsoilmx
330
331day_step = 48 ! default value for day_step
332write(*,*)'Number of time steps per sol?'
333call getin("day_step",day_step)
334write(*,*) " day_step = ",day_step
335
336ecritphy = day_step ! default value for ecritphy, output every time step
337
338ndt = 10 ! default value for ndt
339write(*,*)'Number of sols to run?'
340call getin("ndt",ndt)
341write(*,*) " ndt = ",ndt
342
343dayn = day0 + ndt
344ndt = ndt*day_step
345dttestphys = daysec/day_step
346
347! Imposed surface pressure
348! ------------------------
349psurf = 610. ! Default value for psurf
350write(*,*) 'Surface pressure (Pa)?'
351if (.not. therestart1D) then
352    call getin("psurf",psurf)
353else
354    open(3,file = start1Dname,status = "old",action = "read")
355    read(3,*) header, psurf
356endif
357write(*,*) " psurf = ",psurf
358
359! Reference pressures
360pa = 20.     ! transition pressure (for hybrid coord.)
361preff = 610. ! reference surface pressure
362
363! Aerosol properties
364! ------------------
365tauvis = 0.2 ! default value for tauvis (dust opacity)
366write(*,'("Reference dust opacity at ",f4.0," Pa?")')odpref
367call getin("tauvis",tauvis)
368write(*,*) " tauvis = ",tauvis
369
370! Orbital parameters
371! ------------------
372if (.not. startfiles_1D) then
373    paleomars = .false. ! Default: no water ice reservoir
374    call getin("paleomars",paleomars)
375    if (paleomars) then
376        write(*,*) "paleomars=", paleomars
377        write(*,*) "Orbital parameters from callphys.def"
378        write(*,*) "Enter eccentricity & Lsperi"
379        write(*,*) 'Martian eccentricity (0<e<1)?'
380        call getin('eccentric ',eccentric)
381        write(*,*)"eccentric =",eccentric
382        write(*,*) 'Solar longitude of perihelion (0<Ls<360)?'
383        call getin('Lsperi',Lsperi )
384        write(*,*)"Lsperi=",Lsperi
385        Lsperi = Lsperi*pi/180.0 ! Put it in rad for peri_day
386        periheli = halfaxe*(1 - eccentric)
387        aphelie = halfaxe*(1 + eccentric)
388        call call_dayperi(Lsperi,eccentric,peri_day,year_day)
389        write(*,*) "Corresponding orbital params for GCM"
390        write(*,*) " periheli = ",periheli
391        write(*,*) " aphelie = ",aphelie
392        write(*,*) "date of perihelion (sol)",peri_day
393    else
394        write(*,*) "paleomars=", paleomars
395        write(*,*) "Default present-day orbital parameters"
396        write(*,*) "Unless specified otherwise"
397        write(*,*)'Min. distance Sun-Mars (Mkm)?'
398        call getin("periheli",periheli)
399        write(*,*) " periheli = ",periheli
400
401        write(*,*)'Max. distance Sun-Mars (Mkm)?'
402        call getin("aphelie",aphelie)
403        write(*,*) " aphelie = ",aphelie
404
405        write(*,*)'Day of perihelion?'
406        call getin("periday",peri_day)
407        write(*,*) " periday = ",peri_day
408
409        write(*,*)'Obliquity?'
410        call getin("obliquit",obliquit)
411        write(*,*) " obliquit = ",obliquit
412     endif
413endif !(.not. startfiles_1D )
414
415! Latitude/longitude
416! ------------------
417latitude = 0. ! default value for latitude
418write(*,*)'latitude (in degrees)?'
419call getin("latitude",latitude(1))
420write(*,*) " latitude = ",latitude
421latitude = latitude*pi/180.
422longitude = 0.
423longitude = longitude*pi/180.
424
425! Some initializations (some of which have already been
426! done above!) and loads parameters set in callphys.def
427! and allocates some arrays
428! Mars possible matter with dttestphys in input and include!!!
429! Initializations below should mimick what is done in iniphysiq for 3D GCM
430call init_interface_dyn_phys
431call init_regular_lonlat(1,1,longitude,latitude,(/0.,0./),(/0.,0./))
432call init_geometry(1,longitude,latitude,(/0.,0.,0.,0./),(/0.,0.,0.,0./),cell_area)
433call init_geometry_cell_area_for_outputs(1,cell_area)
434! Ehouarn: init_vertial_layers called later (because disvert not called yet)
435!      call init_vertical_layers(nlayer,preff,scaleheight,
436!     &                      ap,bp,aps,bps,presnivs,pseudoalt)
437call init_dimphy(1,nlayer) ! Initialize dimphy module
438call phys_state_var_init(1,llm,nq,tname,day0,dayn,time,daysec,dttestphys,rad,g,r,cpp,nqperes,nqfils)! MVals: variables isotopes
439call ini_fillgeom(1,latitude,longitude,(/1.0/))
440call conf_phys(1,llm,nq)
441
442! In 1D model physics are called every time step
443! ovverride iphysiq value that has been set by conf_phys
444if (iphysiq /= 1) then
445    write(*,*) "init_testphys1d: setting iphysiq=1"
446    iphysiq = 1
447endif
448
449! Initialize tracers (surface and atmosphere) here:
450write(*,*) "init_testphys1d: initializing tracers"
451if (.not. therestart1D) then
452    call read_profile(nq,nlayer,qsurf(1,:,1),q)
453else
454    do iq = 1,nq
455        read(3,*) header, (qsurf(1,iq,j), j = 1,size(qsurf,3)), (q(1,ilayer,iq), ilayer = 1,nlayer)
456        if (trim(tname(iq)) /= trim(header)) then
457            write(*,*) 'Tracer names not compatible for initialization with "'//trim(start1Dname)//'"!'
458            error stop
459        endif
460    enddo
461endif
462
463! Initialize albedo / soil thermal inertia
464! ----------------------------------------
465if (.not. startfiles_1D) then
466    albedodat(1) = 0.2 ! default value for albedodat
467    write(*,*)'Albedo of bare ground?'
468    call getin("albedo",albedodat(1))
469    write(*,*) " albedo = ",albedodat(1)
470    albedo(1,:,1) = albedodat(1)
471
472    inertiedat(1,1) = 400 ! default value for inertiedat
473    write(*,*)'Soil thermal inertia (SI)?'
474    call getin("inertia",inertiedat(1,1))
475    write(*,*) " inertia = ",inertiedat(1,1)
476
477    ice_depth = -1 ! default value: no ice
478    call getin("subsurface_ice_depth",ice_depth)
479
480    z0(1) = z0_default ! default value for roughness
481    write(*,*) 'Surface roughness length z0 (m)?'
482    call getin("z0",z0(1))
483    write(*,*) " z0 = ",z0(1)
484endif !(.not. startfiles_1D )
485
486! Initialize local slope parameters (only matters if "callslope"
487! is .true. in callphys.def)
488! slope inclination angle (deg) 0: horizontal, 90: vertical
489theta_sl(1) = 0. ! default: no inclination
490call getin("slope_inclination",theta_sl(1))
491! slope orientation (deg)
492! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
493psi_sl(1) = 0. ! default value
494call getin("slope_orientation",psi_sl(1))
495
496! Sub-slopes parameters (assuming no sub-slopes distribution for now).
497def_slope(1) = -90 ! minimum slope angle
498def_slope(2) = 90  ! maximum slope angle
499subslope_dist(1,1) = 1 ! fraction of subslopes in mesh
500
501! For the gravity wave scheme
502! ---------------------------
503zmea(1) = 0.
504zstd(1) = 0.
505zsig(1) = 0.
506zgam(1) = 0.
507zthe(1) = 0.
508
509! For the non-orographic gravity waves scheme
510du_nonoro_gwd(1,:)=0
511dv_nonoro_gwd(1,:)=0
512
513! For the slope wind scheme
514! -------------------------
515hmons(1) = 0.
516write(*,*)'hmons is initialized to ',hmons(1)
517summit(1) = 0.
518write(*,*)'summit is initialized to ',summit(1)
519base(1) = 0.
520
521! Default values initializing the coefficients calculated later
522! -------------------------------------------------------------
523tauscaling(1) = 1.   ! calculated in aeropacity_mod.F
524totcloudfrac(1) = 1. ! calculated in watercloud_mod.F
525
526! Specific initializations for "physiq"
527! -------------------------------------
528! surface geopotential is not used (or useful) since in 1D
529! everything is controled by surface pressure
530phisfi(1) = 0.
531
532! Initialization to take into account prescribed winds
533! ----------------------------------------------------
534ptif = 2.*omeg*sinlat(1)
535
536! Geostrophic wind
537gru = 10. ! default value for gru
538write(*,*)'zonal eastward component of the geostrophic wind (m/s)?'
539call getin("u",gru)
540write(*,*) " u = ",gru
541grv = 0. !default value for grv
542write(*,*)'meridional northward component of the geostrophic wind (m/s)?'
543call getin("v",grv)
544write(*,*) " v = ",grv
545
546! Initialize winds for first time step
547if (.not. therestart1D) then
548    u(:) = gru
549    v(:) = grv
550else
551    read(3,*) header, (u(ilayer), ilayer = 1,nlayer)
552    read(3,*) header, (v(ilayer), ilayer = 1,nlayer)
553endif
554w = 0. ! default: no vertical wind
555
556! Initialize turbulent kinetic energy
557q2 = 0.
558
559! CO2 ice on the surface
560! ----------------------
561! get the index of co2 tracer (not known at this stage)
562igcm_co2 = 0
563do iq = 1,nq
564    if (trim(tname(iq)) == "co2") igcm_co2 = iq
565enddo
566if (igcm_co2 == 0) error stop "init_testphys1d error, missing co2 tracer!"
567
568if (.not. startfiles_1D) then
569    qsurf(1,igcm_co2,1) = 0. ! default value for co2ice
570    write(*,*)'Initial CO2 ice on the surface (kg.m-2)'
571    call getin("co2ice",qsurf(1,igcm_co2,1))
572    write(*,*) " co2ice = ",qsurf(1,igcm_co2,1)
573endif !(.not. startfiles_1D )
574
575! emissivity
576! ----------
577if (.not. startfiles_1D) then
578    emis(:,1) = emissiv
579    if (qsurf(1,igcm_co2,1) == 1.) then
580        emis(:,1) = emisice(1) ! northern hemisphere
581        if (latitude(1) < 0) emis(:,1) = emisice(2) ! southern hemisphere
582    endif
583endif !(.not. startfiles_1D )
584
585! Compute pressures and altitudes of atmospheric levels
586! -----------------------------------------------------
587! Vertical Coordinates
588! """"""""""""""""""""
589hybrid = .true.
590write(*,*)'Hybrid coordinates?'
591call getin("hybrid",hybrid)
592write(*,*) " hybrid = ", hybrid
593
594call disvert_noterre
595! Now that disvert has been called, initialize module vertical_layers_mod
596call init_vertical_layers(nlayer,preff,scaleheight,ap,bp,aps,bps,presnivs,pseudoalt)
597
598plev(:) = ap(:) + psurf*bp(:)
599play(:) = aps(:) + psurf*bps(:)
600zlay(:) = -200.*r*log(play(:)/plev(1))/g
601
602! Initialize temperature profile
603! ------------------------------
604pks = psurf**rcp
605
606! Altitude in km in profile: divide zlay by 1000
607tmp1(0) = 0.
608tmp1(1:) = zlay(:)/1000.
609
610call profile(nlayer + 1,tmp1,tmp2)
611
612if (.not. therestart1D) then
613    tsurf(:,1) = tmp2(0)
614    temp(:) = tmp2(1:)
615else
616    read(3,*) header, (tsurf(1,:), j = 1,size(tsurf,2)), (temp(ilayer), ilayer = 1,nlayer)
617    close(3)
618endif
619
620! Initialize soil properties and temperature
621! ------------------------------------------
622volcapa = 1.e6 ! volumetric heat capacity
623
624if (.not. startfiles_1D) then
625    ! Initialize depths
626    ! -----------------
627    do isoil = 1,nsoil
628        layer(isoil) = 2.e-4*(2.**(isoil - 1)) ! layer depth
629    enddo
630
631    ! Creating the new soil inertia table if there is subsurface ice:
632    if (ice_depth > 0) then
633        iref = 1 ! ice/regolith boundary index
634        if (ice_depth < layer(1)) then
635            inertiedat(1,1) = sqrt(layer(1)/((ice_depth/inertiedat(1,1)**2) + ((layer(1) - ice_depth)/inertieice**2)))
636            inertiedat(1,2:) = inertieice
637        else ! searching for the ice/regolith boundary:
638            do isoil = 1,nsoil
639                if ((ice_depth >= layer(isoil)) .and. (ice_depth < layer(isoil + 1))) then
640                    iref = isoil + 1
641                    exit
642                endif
643            enddo
644            ! We then change the soil inertia table:
645            inertiedat(1,:iref - 1) = inertiedat(1,1)
646            ! We compute the transition in layer(iref)
647            inertiedat(1,iref) = sqrt((layer(iref) - layer(iref - 1))/(((ice_depth - layer(iref - 1))/inertiedat(1,1)**2) + ((layer(iref) - ice_depth)/inertieice**2)))
648            ! Finally, we compute the underlying ice:
649            inertiedat(1,iref + 1:) = inertieice
650        endif ! (ice_depth < layer(1))
651    else ! ice_depth < 0 all is set to surface thermal inertia
652        inertiedat(1,:) = inertiedat(1,1) ! soil thermal inertia
653    endif ! ice_depth > 0
654
655    inertiesoil(1,:,1) = inertiedat(1,:)
656
657    tsoil(:,:,1) = tsurf(1,1) ! soil temperature
658endif !(.not. startfiles_1D)
659
660flux_geo_tmp = 0.
661call getin("flux_geo",flux_geo_tmp)
662flux_geo(:,:) = flux_geo_tmp
663
664! Initialize depths
665! -----------------
666do isoil = 0,nsoil - 1
667    mlayer(isoil) = 2.e-4*(2.**(isoil - 0.5)) ! mid-layer depth
668    layer(isoil + 1) = 2.e-4*(2.**isoil)      ! layer depth
669enddo
670
671! Initialize traceurs
672! -------------------
673if (photochem .or. callthermos) then
674    write(*,*) 'Initializing chemical species'
675    ! flagthermo=0: initialize over all atmospheric layers
676    flagthermo = 0
677    ! check if "h2o_vap" has already been initialized
678    ! (it has been if there is a "profile_h2o_vap" file around)
679    inquire(file = "profile_h2o_vap",exist = there)
680    if (there) then
681        flagh2o = 0 ! 0: do not initialize h2o_vap
682    else
683        flagh2o = 1 ! 1: initialize h2o_vap in inichim_newstart
684    endif
685
686    ! hack to accomodate that inichim_newstart assumes that
687    ! q & psurf arrays are on the dynamics scalar grid
688    allocate(qdyn(2,1,llm,nq),psdyn(2,1))
689    qdyn(1,1,1:llm,1:nq) = q(1,1:llm,1:nq)
690    psdyn(1:2,1) = psurf
691    call inichim_newstart(ngrid,nq,qdyn,qsurf(1,:,1),psdyn,flagh2o,flagthermo)
692    q(1,1:llm,1:nq) = qdyn(1,1,1:llm,1:nq)
693endif
694
695! Check if the surface is a water ice reservoir
696! ---------------------------------------------
697if (.not. startfiles_1D) watercap(1,:) = 0 ! Initialize watercap
698watercaptag(1) = .false. ! Default: no water ice reservoir
699write(*,*)'Water ice cap on ground?'
700call getin("watercaptag",watercaptag)
701write(*,*) " watercaptag = ",watercaptag
702
703! Check if the atmospheric water profile is specified
704! ---------------------------------------------------
705! Adding an option to force atmospheric water values JN
706atm_wat_profile = -1. ! Default: free atm wat profile
707if (water) then
708    write(*,*)'Force atmospheric water vapor profile?'
709    call getin('atm_wat_profile',atm_wat_profile)
710    write(*,*) 'atm_wat_profile = ', atm_wat_profile
711    if (abs(atm_wat_profile + 1.) < 1.e-15) then ! if == -1.
712        write(*,*) 'Free atmospheric water vapor profile'
713        write(*,*) 'Total water is conserved in the column'
714    else if (abs(atm_wat_profile) < 1.e-15) then ! if == 0.
715        write(*,*) 'Dry atmospheric water vapor profile'
716    else if (0. < atm_wat_profile .and. atm_wat_profile <= 1.) then
717        write(*,*) 'Prescribed atmospheric water vapor profile'
718        write(*,*) 'Unless it reaches saturation (maximal value)'
719    else
720        error stop 'Water vapor profile value not correct!'
721    endif
722endif
723
724! Check if the atmospheric water profile relaxation is specified
725! --------------------------------------------------------------
726! Adding an option to relax atmospheric water values JBC
727atm_wat_tau = -1. ! Default: no time relaxation
728if (water) then
729    write(*,*) 'Relax atmospheric water vapor profile?'
730    call getin('atm_wat_tau',atm_wat_tau)
731    write(*,*) 'atm_wat_tau = ', atm_wat_tau
732    if (atm_wat_tau < 0.) then
733        write(*,*) 'Atmospheric water vapor profile is not relaxed.'
734    else
735        if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) then
736            write(*,*) 'Relaxed atmospheric water vapor profile towards ', atm_wat_profile
737            write(*,*) 'Unless it reaches saturation (maximal value)'
738        else
739            write(*,*) 'Reference atmospheric water vapor profile not known!'
740            error stop 'Please, specify atm_wat_profile'
741        endif
742    endif
743endif
744
745END SUBROUTINE init_testphys1d
746
747END MODULE init_testphys1d_mod
Note: See TracBrowser for help on using the repository browser.