source: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F @ 3021

Last change on this file since 3021 was 3021, checked in by jbclement, 16 months ago

Mars PCM 1D:
Addition in 1D of relaxation towards a profile for atmospheric water. It needs the flag "atm_wat_tau" (real):

  • atm_wat_profile < 0. -> no relaxation (default);
  • atm_wat_profile >= 0. -> relaxation towards the value "atm_wat_profile" with relaxation time "atm_wat_tau".
File size: 36.2 KB
Line 
1
2      PROGRAM testphys1d
3! to use  'getin'
4      USE ioipsl_getincom, only: getin
5      use dimphy, only : init_dimphy
6      use mod_grid_phy_lmdz, only : regular_lonlat
7      use infotrac, only: nqtot, tname, nqperes,nqfils
8      use comsoil_h, only: volcapa, layer, mlayer, inertiedat,
9     &                     inertiesoil,nsoilmx,flux_geo
10      use comgeomfi_h, only: sinlat, ini_fillgeom
11      use surfdat_h, only: albedodat, z0_default, emissiv, emisice,
12     &                     albedice, iceradius, dtemisice, z0,
13     &                     zmea, zstd, zsig, zgam, zthe, phisfi,
14     &                     watercaptag, watercap, hmons, summit, base,
15     &                     perenial_co2ice
16      use slope_mod, only: theta_sl, psi_sl
17      use comslope_mod, only: def_slope,subslope_dist,def_slope_mean
18      use phyredem, only: physdem0,physdem1
19      use geometry_mod, only: init_geometry
20      use watersat_mod, only: watersat
21      use tracer_mod, only: igcm_h2o_vap,igcm_h2o_ice,igcm_co2,noms
22      use planete_h, only: year_day, periheli, aphelie, peri_day,
23     &                     obliquit, emin_turb, lmixmin
24      use comcstfi_h, only: pi, rad, omeg, g, mugaz, rcp, r, cpp
25      use time_phylmdz_mod, only: daysec, dtphys, day_step,
26     &                            ecritphy, iphysiq
27      use dimradmars_mod, only: tauvis,totcloudfrac
28      use dust_param_mod, only: tauscaling
29      USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff,sig,
30     &                       presnivs,pseudoalt,scaleheight
31      USE vertical_layers_mod, ONLY: init_vertical_layers
32      USE logic_mod, ONLY: hybrid
33      use physics_distribution_mod, only: init_physics_distribution
34      use regular_lonlat_mod, only: init_regular_lonlat
35      use mod_interface_dyn_phys, only: init_interface_dyn_phys
36      USE phys_state_var_init_mod, ONLY: phys_state_var_init
37      USE physiq_mod, ONLY: physiq
38      USE read_profile_mod, ONLY: read_profile
39      use inichim_newstart_mod, only: inichim_newstart
40      use phyetat0_mod, only: phyetat0
41      USE netcdf, only: NF90_OPEN, NF90_NOERR, NF90_NOWRITE,
42     &                  nf90_strerror
43      use iostart, only: open_startphy,get_var, close_startphy
44      use write_output_mod, only: write_output
45! mostly for  XIOS outputs:
46      use mod_const_mpi, only: init_const_mpi, COMM_LMDZ
47      use parallel_lmdz, only: init_parallel
48
49      IMPLICIT NONE
50
51c=======================================================================
52c   subject:
53c   --------
54c   PROGRAM useful to run physical part of the martian GCM in a 1D column
55c       
56c Can be compiled with a command like (e.g. for 25 layers)
57c  "makegcm -p mars -d 25 testphys1d"
58c It requires the files "testphys1d.def" "callphys.def"
59c   and a 'run.def' file (containing a "INCLUDEDEF=callphys.def" line)
60c      and a file describing the sigma layers (e.g. "z2sig.def")
61c
62c   author: Frederic Hourdin, R.Fournier,F.Forget
63c   -------
64c   
65c   update: 12/06/2003 including chemistry (S. Lebonnois)
66c                            and water ice (F. Montmessin)
67c
68c=======================================================================
69
70      include "dimensions.h"
71      integer, parameter :: ngrid = 1 !(2+(jjm-1)*iim - 1/jjm)
72      integer, parameter :: nlayer = llm
73!#include "dimradmars.h"
74!#include "comgeomfi.h"
75!#include "surfdat.h"
76!#include "slope.h"
77!#include "comsoil.h"
78!#include "comdiurn.h"
79      include "callkeys.h"
80!#include "comsaison.h"
81!#include "control.h"
82      include "netcdf.inc"
83!#include "advtrac.h"
84
85c --------------------------------------------------------------
86c  Declarations
87c --------------------------------------------------------------
88c
89      INTEGER unitstart      ! unite d'ecriture de "startfi"
90      INTEGER nlevel,nsoil,ndt
91      INTEGER ilayer,ilevel,isoil,idt,iq
92      LOGICAl firstcall,lastcall
93      LOGICAL write_prof
94c
95      real,parameter :: odpref=610. ! DOD reference pressure (Pa)
96c
97      INTEGER day0,dayn          ! date initial (sol ; =0 a Ls=0) and final
98      REAL day           ! date durant le run
99      REAL time             ! time (0<time<1 ; time=0.5 a midi)
100      REAL play(nlayer)   ! Pressure at the middle of the layers (Pa)
101      REAL plev(nlayer+1) ! intermediate pressure levels (pa)
102      REAL psurf,tsurf(1)     
103      REAL u(nlayer),v(nlayer)  ! zonal, meridional wind
104      REAL gru,grv   ! prescribed "geostrophic" background wind
105      REAL temp(nlayer)   ! temperature at the middle of the layers
106      REAL,ALLOCATABLE :: q(:,:) ! tracer mixing ratio (e.g. kg/kg)
107      REAL,ALLOCATABLE :: qsurf(:) ! tracer surface budget (e.g. kg.m-2)
108      REAL tsoil(nsoilmx)   ! subsurface soik temperature (K)
109      REAL emis(1)          ! surface layer
110      REAL albedo(1,1)      ! surface albedo
111      REAL :: wstar(1)=0.    ! Thermals vertical velocity
112      REAL q2(nlayer+1)   ! Turbulent Kinetic Energy
113      REAL zlay(nlayer)   ! altitude estimee dans les couches (km)
114
115c    Physical and dynamical tandencies (e.g.  m.s-2, K/s, Pa/s)
116      REAL du(nlayer),dv(nlayer),dtemp(nlayer)
117      REAL dudyn(nlayer),dvdyn(nlayer),dtempdyn(nlayer)
118      REAL dpsurf(1)   
119      REAL,ALLOCATABLE :: dq(:,:)
120      REAL,ALLOCATABLE :: dqdyn(:,:)
121
122c   Various intermediate variables
123      INTEGER flagthermo,flagh2o
124      REAL zls
125      REAL phi(nlayer),h(nlayer),s(nlayer)
126      REAL pks, ptif, w(nlayer)
127      REAL qtotinit,qtot
128      INTEGER ierr, aslun
129      REAL tmp1(0:nlayer),tmp2(0:nlayer)
130      integer :: nq=1 ! number of tracers
131      real :: latitude(1), longitude(1), cell_area(1)
132      ! dummy variables along "dynamics scalar grid"
133      real,allocatable :: qdyn(:,:,:,:),psdyn(:,:)
134
135      character*2 str2
136      character (len=7) :: str7
137      character(len=44) :: txt
138
139c   New flag to compute paleo orbital configurations + few variables JN
140      REAL halfaxe, excentric, Lsperi
141      Logical paleomars
142c   JN : Force atmospheric water profiles
143      REAL atm_wat_profile
144      REAL atm_wat_tau
145      REAL,ALLOCATABLE :: zqsat(:) ! useful for (atm_wat_profile=2)
146
147
148c   MVals: isotopes as in the dynamics (CRisi)
149      INTEGER :: ifils,ipere,generation
150      CHARACTER(len=30), ALLOCATABLE, DIMENSION(:) :: tnom_transp ! transporting fluid short name
151      CHARACTER(len=80) :: line ! to store a line of text     
152      INTEGER ierr0
153      LOGICAL :: continu
154
155      logical :: present
156
157c   LL: Subsurface geothermal flux
158      real :: flux_geo_tmp
159
160c   RV: Start from a startfi and not run.def
161      logical :: startfile_1D
162      REAL tab_cntrl(100)
163      LOGICAL :: found
164      REAL albedo_read(1,2,1)      ! surface albedo
165
166c   LL: Possibility to add subsurface ice
167      REAL :: ice_depth            ! depth of the ice table, ice_depth < 0. means no ice
168      REAL :: inertieice = 2100.   ! ice thermal inertia
169      integer :: iref
170 
171c=======================================================================
172
173c=======================================================================
174c INITIALISATION
175c=======================================================================
176#ifdef CPP_XIOS
177      CALL init_const_mpi
178      CALL init_parallel
179#endif
180
181! initialize "serial/parallel" related stuff
182!      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
183!      CALL init_phys_lmdz(1,1,llm,1,(/1/))
184!      call initcomgeomphy
185
186
187c ------------------------------------------------------
188c  Loading run parameters from "run.def" file
189c ------------------------------------------------------
190
191
192! check if 'run.def' file is around (otherwise reading parameters
193! from callphys.def via getin() routine won't work.
194      open(99,file='run.def',status='old',form='formatted',
195     &     iostat=ierr)
196      if (ierr.ne.0) then
197        write(*,*) 'Cannot find required file "run.def"'
198        write(*,*) '  (which should contain some input parameters'
199        write(*,*) '   along with the following line:'
200        write(*,*) '   INCLUDEDEF=callphys.def'
201        write(*,*) '   )'
202        write(*,*) ' ... might as well stop here ...'
203        stop
204      else
205        close(99)
206      endif
207
208      write(*,*)'Do you want to start with a startfi and write
209     & a restartfi?'
210      startfile_1D=.false.
211      call getin("startfile_1D",startfile_1D)
212      write(*,*)" startfile_1D = ", startfile_1D
213
214c ------------------------------------------------------
215c  Prescribed constants to be set here
216c ------------------------------------------------------
217
218c      if(.not. startfile_1D ) then
219
220      pi=2.E+0*asin(1.E+0)
221
222c     Mars planetary constants
223c     ----------------------------
224      rad=3397200.               ! mars radius (m)  ~3397200 m
225      daysec=88775.              ! length of a sol (s)  ~88775 s
226      omeg=4.*asin(1.)/(daysec)  ! rotation rate (rad.s-1)
227      g=3.72                     ! gravity (m.s-2) ~3.72 
228      mugaz=43.49                ! atmosphere mola mass (g.mol-1) ~43.49
229      rcp=.256793                ! = r/cp  ~0.256793
230      r= 8.314511E+0 *1000.E+0/mugaz
231      cpp= r/rcp
232      year_day = 669             ! lenght of year (sols) ~668.6
233      periheli = 206.66          ! minimum sun-mars distance (Mkm) ~206.66
234      aphelie = 249.22           ! maximum sun-mars distance (Mkm) ~249.22
235      halfaxe = 227.94           ! demi-grand axe de l'ellipse
236      peri_day =  485.           ! perihelion date (sols since N. Spring)
237      obliquit = 25.2            ! Obliquity (deg) ~25.2         
238      excentric = 0.0934         ! Eccentricity (0.0934)         
239 
240c     Planetary Boundary Layer and Turbulence parameters
241c     --------------------------------------------------
242      z0_default =  1.e-2        ! surface roughness (m) ~0.01
243      emin_turb = 1.e-6          ! minimal turbulent energy ~1.e-8
244      lmixmin = 30               ! mixing length ~100
245 
246c     cap properties and surface emissivities
247c     ----------------------------------------------------
248      emissiv= 0.95              ! Bare ground emissivity ~.95
249      emisice(1)=0.95            ! Northern cap emissivity
250      emisice(2)=0.95            ! Southern cap emisssivity
251      albedice(1)=0.5            ! Northern cap albedo
252      albedice(2)=0.5            ! Southern cap albedo
253      iceradius(1) = 100.e-6     ! mean scat radius of CO2 snow (north)
254      iceradius(2) = 100.e-6     ! mean scat radius of CO2 snow (south)
255      dtemisice(1) = 2.          ! time scale for snow metamorphism (north)
256      dtemisice(2) = 2.          ! time scale for snow metamorphism (south)
257
258computeTcondc      endif ! .not. startfile_1D
259
260c     mesh surface (not a very usefull quantity in 1D)
261c     ----------------------------------------------------
262      cell_area(1)=1.E+0
263
264
265! check if there is a 'traceur.def' file
266! and process it.
267      ! load tracer names from file 'traceur.def'
268        open(90,file='traceur.def',status='old',form='formatted',
269     &       iostat=ierr)
270        if (ierr.ne.0) then
271          write(*,*) 'Cannot find required file "traceur.def"'
272          write(*,*) ' If you want to run with tracers, I need it'
273          write(*,*) ' ... might as well stop here ...'
274          stop
275        else
276          write(*,*) "testphys1d: Reading file traceur.def"
277          ! read number of tracers:
278          read(90,*,iostat=ierr) nq
279          nqtot=nq ! set value of nqtot (in infotrac module) as nq
280          if (ierr.ne.0) then
281            write(*,*) "testphys1d: error reading number of tracers"
282            write(*,*) "   (first line of traceur.def) "
283            stop
284          endif
285          if (nq<1) then
286            write(*,*) "testphys1d: error number of tracers"
287            write(*,*) "is nq=",nq," but must be >=1!"
288            stop
289          endif
290        endif
291        ! allocate arrays:
292        allocate(tname(nq))
293        allocate(q(nlayer,nq))
294        allocate(zqsat(nlayer))
295        allocate(qsurf(nq))
296        allocate(dq(nlayer,nq))
297        allocate(dqdyn(nlayer,nq))
298        allocate(tnom_transp(nq))
299       
300        ! read tracer names from file traceur.def
301        do iq=1,nq
302          read(90,'(80a)',iostat=ierr) line ! store the line from traceur.def
303          if (ierr.ne.0) then
304            write(*,*) 'testphys1d: error reading tracer names...'
305            stop
306          endif
307          ! if format is tnom_0, tnom_transp (isotopes)
308          read(line,*,iostat=ierr0) tname(iq),tnom_transp(iq)
309          if (ierr0.ne.0) then
310            read(line,*) tname(iq)
311            tnom_transp(iq)='air'
312          endif
313
314        enddo
315        close(90)
316
317       ! Isotopes: as in the 3D case we have to determine father/son relations for isotopes and carrying fluid
318       ALLOCATE(nqfils(nqtot))
319       nqperes=0
320       nqfils(:)=0 
321       DO iq=1,nqtot
322       if (tnom_transp(iq) == 'air') then
323         ! ceci est un traceur père
324         WRITE(*,*) 'Le traceur',iq,', appele ',
325     &          trim(tname(iq)),', est un pere'
326         nqperes=nqperes+1
327       else !if (tnom_transp(iq) == 'air') then
328         ! ceci est un fils. Qui est son père?
329         WRITE(*,*) 'Le traceur',iq,', appele ',
330     &                trim(tname(iq)),', est un fils'
331         continu=.true.
332         ipere=1
333         do while (continu)           
334           if (tnom_transp(iq) .eq. tname(ipere)) then
335             ! Son père est ipere
336             WRITE(*,*) 'Le traceur',iq,'appele ',
337     &   trim(tname(iq)),' est le fils de ',
338     &   ipere,'appele ',trim(tname(ipere))
339             nqfils(ipere)=nqfils(ipere)+1         
340             continu=.false.
341           else !if (tnom_transp(iq) == tnom_0(ipere)) then
342             ipere=ipere+1
343             if (ipere.gt.nqtot) then
344                 WRITE(*,*) 'Le traceur',iq,'appele ',
345     &           trim(tname(iq)),', est orpelin.'
346                 CALL abort_gcm('infotrac_init',
347     &                  'Un traceur est orphelin',1)
348             endif !if (ipere.gt.nqtot) then
349           endif !if (tnom_transp(iq) == tnom_0(ipere)) then
350         enddo !do while (continu)
351       endif !if (tnom_transp(iq) == 'air') then
352       enddo !DO iq=1,nqtot
353       WRITE(*,*) 'nqperes=',nqperes   
354       WRITE(*,*) 'nqfils=',nqfils
355
356        ! initialize tracers here:
357        write(*,*) "testphys1d: initializing tracers"
358        call read_profile(nq, nlayer, qsurf, q)
359
360      call init_physics_distribution(regular_lonlat,4,
361     &                               1,1,1,nlayer,COMM_LMDZ)
362
363      if(.not. startfile_1D ) then
364
365c  Date and local time at beginning of run
366c  ---------------------------------------
367c    Date (in sols since spring solstice) at beginning of run
368      day0 = 0 ! default value for day0
369      write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?'
370      call getin("day0",day0)
371      day=float(day0)
372      write(*,*) " day0 = ",day0
373c  Local time at beginning of run
374      time=0 ! default value for time
375      write(*,*)'Initial local time (in hours, between 0 and 24)?'
376      call getin("time",time)
377      write(*,*)" time = ",time
378      time=time/24.E+0 ! convert time (hours) to fraction of sol
379
380      else
381      call open_startphy("startfi.nc")
382      call get_var("controle",tab_cntrl,found)
383       if (.not.found) then
384         call abort_physic("open_startphy",
385     &        "tabfi: Failed reading <controle> array",1)
386       else
387         write(*,*)'tabfi: tab_cntrl',tab_cntrl
388       endif
389       day0 = tab_cntrl(3)
390       day=float(day0)
391
392       call get_var("Time",time,found)
393
394       call close_startphy
395
396      endif !startfile_1D
397
398c  Discretization (Definition of grid and time steps)
399c  --------------
400c
401      nlevel=nlayer+1
402      nsoil=nsoilmx
403
404      day_step=48 ! default value for day_step
405      PRINT *,'Number of time steps per sol ?'
406      call getin("day_step",day_step)
407      write(*,*) " day_step = ",day_step
408
409      ecritphy=day_step ! default value for ecritphy, output every time step
410
411      ndt=10 ! default value for ndt
412      PRINT *,'Number of sols to run ?'
413      call getin("ndt",ndt)
414      write(*,*) " ndt = ",ndt
415
416      dayn=day0+ndt
417      ndt=ndt*day_step     
418      dtphys=daysec/day_step 
419
420c Imposed surface pressure
421c ------------------------------------
422c
423
424      psurf=610. ! default value for psurf
425      PRINT *,'Surface pressure (Pa) ?'
426      call getin("psurf",psurf)
427      write(*,*) " psurf = ",psurf
428
429c Reference pressures
430      pa=20.   ! transition pressure (for hybrid coord.)
431      preff=610.      ! reference surface pressure
432 
433c Aerosol properties
434c --------------------------------
435      tauvis=0.2 ! default value for tauvis (dust opacity)
436      write(*,'("Reference dust opacity at ",f4.0," Pa ?")')odpref
437      call getin("tauvis",tauvis)
438      write(*,*) " tauvis = ",tauvis
439
440c Orbital parameters
441c ------------------
442
443      if(.not. startfile_1D ) then
444
445      paleomars=.false. ! Default: no water ice reservoir
446      call getin("paleomars",paleomars)
447      if (paleomars.eqv..true.) then
448        write(*,*) "paleomars=", paleomars
449        write(*,*) "Orbital parameters from callphys.def"
450        write(*,*) "Enter eccentricity & Lsperi"
451        print *, 'Martian eccentricity (0<e<1) ?'
452        call getin('excentric ',excentric)
453        write(*,*)"excentric =",excentric
454        print *, 'Solar longitude of perihelion (0<Ls<360) ?'
455        call getin('Lsperi',Lsperi )
456        write(*,*)"Lsperi=",Lsperi
457        Lsperi = Lsperi*pi/180.0 ! Put it in rad for peri_day
458        periheli = halfaxe*(1-excentric)
459        aphelie = halfaxe*(1+excentric)
460        call call_dayperi(Lsperi,excentric,peri_day,year_day)
461        write(*,*) "Corresponding orbital params for GCM"
462        write(*,*) " periheli = ",periheli
463        write(*,*) " aphelie = ",aphelie
464        write(*,*) "date of perihelion (sol)",peri_day
465      else
466        write(*,*) "paleomars=", paleomars
467        write(*,*) "Default present-day orbital parameters"
468        write(*,*) "Unless specified otherwise"
469        print *,'Min. distance Sun-Mars (Mkm)?'
470        call getin("periheli",periheli)
471        write(*,*) " periheli = ",periheli
472
473        print *,'Max. distance Sun-Mars (Mkm)?'
474        call getin("aphelie",aphelie)
475        write(*,*) " aphelie = ",aphelie
476
477        print *,'Day of perihelion?'
478        call getin("periday",peri_day)
479        write(*,*) " periday = ",peri_day
480
481        print *,'Obliquity?'
482        call getin("obliquit",obliquit)
483        write(*,*) " obliquit = ",obliquit
484      end if
485
486      endif !(.not. startfile_1D )
487 
488c  latitude/longitude
489c  ------------------
490      latitude(1)=0 ! default value for latitude
491      PRINT *,'latitude (in degrees) ?'
492      call getin("latitude",latitude(1))
493      write(*,*) " latitude = ",latitude
494      latitude=latitude*pi/180.E+0
495      longitude=0.E+0
496      longitude=longitude*pi/180.E+0
497
498!  some initializations (some of which have already been
499!  done above!) and loads parameters set in callphys.def
500!  and allocates some arrays
501!Mars possible matter with dtphys in input and include!!!
502! Initializations below should mimick what is done in iniphysiq for 3D GCM
503      call init_interface_dyn_phys
504      call init_regular_lonlat(1,1,longitude,latitude,
505     &                   (/0.,0./),(/0.,0./))
506      call init_geometry(1,longitude,latitude,
507     &                   (/0.,0.,0.,0./),(/0.,0.,0.,0./),
508     &                   cell_area)
509! Ehouarn: init_vertial_layers called later (because disvert not called yet)
510!      call init_vertical_layers(nlayer,preff,scaleheight,
511!     &                      ap,bp,aps,bps,presnivs,pseudoalt)
512      call init_dimphy(1,nlayer) ! Initialize dimphy module
513      call phys_state_var_init(1,llm,nq,tname,
514     .          day0,dayn,time,
515     .          daysec,dtphys,
516     .          rad,g,r,cpp,
517     .          nqperes,nqfils)! MVals: variables isotopes
518      call ini_fillgeom(1,latitude,longitude,(/1.0/))
519      call conf_phys(1,llm,nq)
520
521      ! in 1D model physics are called every time step
522      ! ovverride iphysiq value that has been set by conf_phys
523      if (iphysiq/=1) then
524        write(*,*) "testphys1d: setting iphysiq=1"
525        iphysiq=1
526      endif
527
528c  Initialize albedo / soil thermal inertia
529c  ----------------------------------------
530c
531
532      if(.not. startfile_1D ) then
533
534      albedodat(1)=0.2 ! default value for albedodat
535      PRINT *,'Albedo of bare ground ?'
536      call getin("albedo",albedodat(1))
537      write(*,*) " albedo = ",albedodat(1)
538      albedo(1,1)=albedodat(1)
539
540      inertiedat(1,1)=400 ! default value for inertiedat
541      PRINT *,'Soil thermal inertia (SI) ?'
542      call getin("inertia",inertiedat(1,1))
543      write(*,*) " inertia = ",inertiedat(1,1)
544
545      ice_depth = -1 ! default value: no ice
546      CALL getin("subsurface_ice_depth",ice_depth)
547
548      z0(1)=z0_default ! default value for roughness
549      write(*,*) 'Surface roughness length z0 (m)?'
550      call getin("z0",z0(1))
551      write(*,*) " z0 = ",z0(1)
552
553      endif !(.not. startfile_1D )
554
555! Initialize local slope parameters (only matters if "callslope"
556! is .true. in callphys.def)
557      ! slope inclination angle (deg) 0: horizontal, 90: vertical
558      theta_sl(1)=0.0 ! default: no inclination
559      call getin("slope_inclination",theta_sl(1))
560      ! slope orientation (deg)
561      ! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
562      psi_sl(1)=0.0 ! default value
563      call getin("slope_orientation",psi_sl(1))
564     
565      ! sub-slopes parameters (assuming no sub-slopes distribution for now).
566      def_slope(1)=-90 ! minimum slope angle
567      def_slope(2)=90 ! maximum slope angle
568      subslope_dist(1,1)=1 ! fraction of subslopes in mesh
569c
570c  for the gravity wave scheme
571c  ---------------------------------
572c
573      zmea(1)=0.E+0
574      zstd(1)=0.E+0
575      zsig(1)=0.E+0
576      zgam(1)=0.E+0
577      zthe(1)=0.E+0
578c
579c  for the slope wind scheme
580c  ---------------------------------
581
582      hmons(1)=0.E+0
583      PRINT *,'hmons is initialized to ',hmons(1)
584      summit(1)=0.E+0
585      PRINT *,'summit is initialized to ',summit(1)
586      base(1)=0.E+0
587c
588c  Default values initializing the coefficients calculated later
589c  ---------------------------------
590
591      tauscaling(1)=1. ! calculated in aeropacity_mod.F
592      totcloudfrac(1)=1. ! calculated in watercloud_mod.F     
593
594c   Specific initializations for "physiq"
595c   -------------------------------------
596c   surface geopotential is not used (or useful) since in 1D
597c   everything is controled by surface pressure
598      phisfi(1)=0.E+0
599
600c   Initialization to take into account prescribed winds
601c   ------------------------------------------------------
602      ptif=2.E+0*omeg*sinlat(1)
603 
604c    geostrophic wind
605      gru=10. ! default value for gru
606      PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
607      call getin("u",gru)
608      write(*,*) " u = ",gru
609      grv=0. !default value for grv
610      PRINT *,'meridional northward component of the geostrophic',
611     &' wind (m/s) ?'
612      call getin("v",grv)
613      write(*,*) " v = ",grv
614
615c     Initialize winds  for first time step
616      DO ilayer=1,nlayer
617         u(ilayer)=gru
618         v(ilayer)=grv
619         w(ilayer)=0 ! default: no vertical wind
620      ENDDO
621
622c     Initialize turbulente kinetic energy
623      DO ilevel=1,nlevel
624         q2(ilevel)=0.E+0
625      ENDDO
626
627c  CO2 ice on the surface
628c  -------------------
629      ! get the index of co2 tracer (not known at this stage)
630      igcm_co2=0
631      do iq=1,nq
632        if (trim(tname(iq))=="co2") then
633          igcm_co2=iq
634        endif
635      enddo
636      if (igcm_co2==0) then
637        write(*,*) "testphys1d error, missing co2 tracer!"
638        stop
639      endif
640
641      if(.not. startfile_1D ) then
642      qsurf(igcm_co2)=0.E+0 ! default value for co2ice
643      PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
644      call getin("co2ice",qsurf(igcm_co2))
645      write(*,*) " co2ice = ",qsurf(igcm_co2)
646      endif !(.not. startfile_1D )
647
648c
649c  emissivity
650c  ----------
651      if(.not. startfile_1D ) then
652      emis=emissiv
653      IF (qsurf(igcm_co2).eq.1.E+0) THEN
654         emis=emisice(1) ! northern hemisphere
655         IF(latitude(1).LT.0) emis=emisice(2) ! southern hemisphere
656      ENDIF
657      endif !(.not. startfile_1D )
658 
659
660c  Compute pressures and altitudes of atmospheric levels
661c  ----------------------------------------------------------------
662
663c    Vertical Coordinates
664c    """"""""""""""""""""
665      hybrid=.true.
666      PRINT *,'Hybrid coordinates ?'
667      call getin("hybrid",hybrid)
668      write(*,*) " hybrid = ", hybrid
669
670      CALL  disvert_noterre
671      ! now that disvert has been called, initialize module vertical_layers_mod
672      call init_vertical_layers(nlayer,preff,scaleheight,
673     &                      ap,bp,aps,bps,presnivs,pseudoalt)
674
675      DO ilevel=1,nlevel
676        plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
677      ENDDO
678
679      DO ilayer=1,nlayer
680        play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
681      ENDDO
682
683      DO ilayer=1,nlayer
684        zlay(ilayer)=-200.E+0 *r*log(play(ilayer)/plev(1))
685     &   /g
686      ENDDO
687
688
689c  Initialize temperature profile
690c  --------------------------------------
691      pks=psurf**rcp
692
693c altitude in km in profile: divide zlay by 1000
694      tmp1(0)=0.E+0
695      DO ilayer=1,nlayer
696        tmp1(ilayer)=zlay(ilayer)/1000.E+0
697      ENDDO
698
699      call profile(nlayer+1,tmp1,tmp2)
700
701      tsurf=tmp2(0)
702      DO ilayer=1,nlayer
703        temp(ilayer)=tmp2(ilayer)
704      ENDDO
705     
706
707! Initialize soil properties and temperature
708! ------------------------------------------
709      volcapa=1.e6 ! volumetric heat capacity
710
711      if(.not. startfile_1D ) then
712
713! Initialize depths
714! -----------------
715       do isoil=1,nsoil
716         layer(isoil)=2.e-4*(2.**(isoil-1)) ! layer depth
717       enddo
718
719! Creating the new soil inertia table if there is subsurface ice :
720       IF (ice_depth.gt.0) THEN
721         iref = 1 ! ice/regolith boundary index
722           IF (ice_depth.lt.layer(1)) THEN
723             inertiedat(1,1) = sqrt( layer(1) /
724     &        ( (ice_depth/inertiedat(1,1)**2) +
725     &        ((layer(1)-ice_depth)/inertieice**2) ) )
726             DO isoil=1,nsoil
727              inertiedat(1,isoil) = inertieice
728             ENDDO   
729           ELSE ! searching for the ice/regolith boundary :
730           DO isoil=1,nsoil
731            IF ((ice_depth.ge.layer(isoil)).and.
732     &         (ice_depth.lt.layer(isoil+1))) THEN
733                  iref=isoil+1
734                  EXIT
735            ENDIF
736           ENDDO
737!         We then change the soil inertia table :
738           DO isoil=1,iref-1
739            inertiedat(1,isoil) = inertiedat(1,1)
740           ENDDO
741!         We compute the transition in layer(iref)
742           inertiedat(1,iref) = sqrt( (layer(iref)-layer(iref-1)) /
743     &          ( ((ice_depth-layer(iref-1))/inertiedat(1,1)**2) +
744     &          ((layer(iref)-ice_depth)/inertieice**2) ) )
745!         Finally, we compute the underlying ice :
746           DO isoil=iref+1,nsoil
747            inertiedat(1,isoil) = inertieice
748           ENDDO
749         ENDIF ! (ice_depth.lt.layer(1))     
750       ELSE ! ice_depth <0 all is set to surface thermal inertia
751         DO isoil=1,nsoil
752           inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
753         ENDDO
754       ENDIF ! ice_depth.gt.0
755
756       inertiesoil(1,:,1) = inertiedat(1,:)
757
758       DO isoil = 1,nsoil
759         tsoil(isoil)=tsurf(1)  ! soil temperature
760       ENDDO
761
762      endif !(.not. startfile_1D )
763
764      flux_geo_tmp=0.
765      call getin("flux_geo",flux_geo_tmp)
766      flux_geo(:,:) = flux_geo_tmp
767
768! Initialize depths
769! -----------------
770      do isoil=0,nsoil-1
771        mlayer(isoil)=2.e-4*(2.**(isoil-0.5)) ! mid-layer depth
772      enddo
773      do isoil=1,nsoil
774        layer(isoil)=2.e-4*(2.**(isoil-1)) ! layer depth
775      enddo
776
777c    Initialize traceurs
778c    ---------------------------
779
780      if (photochem.or.callthermos) then
781         write(*,*) 'Initializing chemical species'
782         ! flagthermo=0: initialize over all atmospheric layers
783         flagthermo=0
784         ! check if "h2o_vap" has already been initialized
785         ! (it has been if there is a "profile_h2o_vap" file around)
786         inquire(file="profile_h2o_vap",exist=present)
787         if (present) then
788           flagh2o=0 ! 0: do not initialize h2o_vap
789         else
790           flagh2o=1 ! 1: initialize h2o_vap in inichim_newstart
791         endif
792         
793         ! hack to accomodate that inichim_newstart assumes that
794         ! q & psurf arrays are on the dynamics scalar grid
795         allocate(qdyn(2,1,llm,nq),psdyn(2,1))
796         qdyn(1,1,1:llm,1:nq)=q(1:llm,1:nq)
797         psdyn(1:2,1)=psurf
798         call inichim_newstart(ngrid,nq,qdyn,qsurf,psdyn,
799     $                         flagh2o,flagthermo)
800         q(1:llm,1:nq)=qdyn(1,1,1:llm,1:nq)
801      endif
802
803c Check if the surface is a water ice reservoir
804c --------------------------------------------------
805      if(.not. startfile_1D ) then
806      watercap(1,:)=0 ! Initialize watercap
807      endif !(.not. startfile_1D )
808      watercaptag(1)=.false. ! Default: no water ice reservoir
809      print *,'Water ice cap on ground ?'
810      call getin("watercaptag",watercaptag)
811      write(*,*) " watercaptag = ",watercaptag
812     
813c Check if the atmospheric water profile is specified
814c ---------------------------------------------------
815      ! Adding an option to force atmospheric water values JN
816      atm_wat_profile = -1. ! Default: free atm wat profile
817      print *,'Force atmospheric water vapor profile?'
818      call getin("atm_wat_profile",atm_wat_profile)
819      write(*,*) "atm_wat_profile = ", atm_wat_profile
820      if (abs(atm_wat_profile + 1.) < 1.e-15) then ! if == -1.
821        write(*,*) "Free atmospheric water vapor profile"
822        write(*,*) "Total water is conserved in the column"
823      else if (abs(atm_wat_profile) < 1.e-15) then ! if == 0.
824        write(*,*) "Dry atmospheric water vapor profile"
825      else if (0. < atm_wat_profile .and. atm_wat_profile <= 1.) then
826        write(*,*) "Prescribed atmospheric water vapor profile"
827        write(*,*) "Unless it reaches saturation (maximal value)"
828      else
829        write(*,*) 'Water vapor profile value not correct!'
830        stop
831      endif
832
833! Check if the atmospheric water profile relaxation is specified
834! --------------------------------------------------------------
835      ! Adding an option to relax atmospheric water values JBC
836      atm_wat_tau = -1. ! Default: no time relaxation
837      print*, 'Relax atmospheric water vapor profile?'
838      call getin("atm_wat_tau",atm_wat_tau)
839      write(*,*) "atm_wat_tau = ", atm_wat_tau
840      if (atm_wat_tau < 0.) then
841        write(*,*) "Atmospheric water vapor profile is not relaxed"
842      else
843        if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) then
844          write(*,*) "Relaxed atmospheric water vapor profile towards ",
845     &    atm_wat_profile
846          write(*,*) "Unless it reaches saturation (maximal value)"
847        else
848       write(*,*) 'Reference atmospheric water vapor profile not known!'
849       write(*,*) 'Please, specify atm_wat_profile'
850       stop
851        endif
852      endif
853
854c  Write a "startfi" file
855c  --------------------
856c  This file will be read during the first call to "physiq".
857c  It is needed to transfert physics variables to "physiq"...
858
859      if(.not. startfile_1D ) then
860
861      call physdem0("startfi.nc",longitude,latitude,nsoilmx,ngrid,
862     &              llm,nq,dtphys,float(day0),0.,cell_area,
863     &              albedodat,inertiedat,def_slope,subslope_dist)
864      call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,
865     &              dtphys,time,
866     &              tsurf,tsoil,inertiesoil,albedo,emis,
867     &              q2,qsurf,tauscaling,
868     &              totcloudfrac,wstar,watercap,perenial_co2ice)
869      endif !(.not. startfile_1D )
870
871c=======================================================================
872c  1D MODEL TIME STEPPING LOOP
873c=======================================================================
874c
875      firstcall=.true.
876      lastcall=.false.
877      DO idt=1,ndt
878        IF (idt.eq.ndt) lastcall=.true.
879c        IF (idt.eq.ndt-day_step-1) then       !test
880c         lastcall=.true.
881c         call solarlong(day*1.0,zls)
882c         write(103,*) 'Ls=',zls*180./pi
883c         write(103,*) 'Lat=', latitude(1)*180./pi
884c         write(103,*) 'Tau=', tauvis/odpref*psurf
885c         write(103,*) 'RunEnd - Atmos. Temp. File'
886c         write(103,*) 'RunEnd - Atmos. Temp. File'
887c         write(104,*) 'Ls=',zls*180./pi
888c         write(104,*) 'Lat=', latitude(1)
889c         write(104,*) 'Tau=', tauvis/odpref*psurf
890c         write(104,*) 'RunEnd - Atmos. Temp. File'
891c        ENDIF
892
893c     compute geopotential
894c     ~~~~~~~~~~~~~~~~~~~~~
895      DO ilayer=1,nlayer
896        s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
897        h(ilayer)=cpp*temp(ilayer)/(pks*s(ilayer))
898      ENDDO
899      phi(1)=pks*h(1)*(1.E+0-s(1))
900      DO ilayer=2,nlayer
901         phi(ilayer)=phi(ilayer-1)+
902     &               pks*(h(ilayer-1)+h(ilayer))*.5E+0
903     &                  *(s(ilayer-1)-s(ilayer))
904
905      ENDDO
906
907!       Modify atmospheric water if asked
908!       Added "atm_wat_profile" flag (JN + JBC)
909!       ---------------------------------------
910        call watersat(nlayer,temp,play,zqsat)
911        if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) then
912        ! If atmospheric water is monitored
913          if (atm_wat_tau < 0.) then ! Prescribed atm_wat_profile
914            ! Wet if >0, Dry if =0
915            q(:,igcm_h2o_vap) = min(zqsat(:),atm_wat_profile*g/psurf)
916            q(:,igcm_h2o_ice) = 0. ! reset h2o ice
917          else ! Relaxation towards the value atm_wat_profile with relaxation time atm_wat_tau
918        q(:,igcm_h2o_vap) = atm_wat_profile*g/psurf + (q(:,igcm_h2o_vap)
919     &              - atm_wat_profile*g/psurf)*dexp(-dtphys/atm_wat_tau)
920            q(:,igcm_h2o_vap) = min(zqsat(:),q(:,igcm_h2o_vap))
921            q(:,igcm_h2o_ice) = 0. ! reset h2o ice
922          endif
923        endif
924
925!      write(*,*) "testphys1d avant q", q(1,:)
926c       call physics
927c       --------------------
928      CALL physiq (1,llm,nq,
929     ,     firstcall,lastcall,
930     ,     day,time,dtphys,
931     ,     plev,play,phi,
932     ,     u, v,temp, q, 
933     ,     w,
934C - outputs
935     s     du, dv, dtemp, dq,dpsurf)
936!      write(*,*) "testphys1d apres q", q(1,:)
937
938c       wind increment : specific for 1D
939c       --------------------------------
940c       The physics compute the tendencies on u and v,
941c       here we just add Coriolos effect
942c
943c       DO ilayer=1,nlayer
944c          du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
945c          dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
946c       ENDDO
947
948c       For some tests : No coriolis force at equator
949c       if(latitude(1).eq.0.) then
950          DO ilayer=1,nlayer
951             du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
952             dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
953          ENDDO
954c       end if
955
956c       Compute time for next time step
957c       ---------------------------------------
958        firstcall=.false.
959        time=time+dtphys/daysec
960        IF (time.gt.1.E+0) then
961            time=time-1.E+0
962            day=day+1
963        ENDIF
964
965c       compute winds and temperature for next time step
966c       ----------------------------------------------------------
967
968        DO ilayer=1,nlayer
969           u(ilayer)=u(ilayer)+dtphys*du(ilayer)
970           v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
971           temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
972        ENDDO
973
974c       compute pressure for next time step
975c       ----------------------------------------------------------
976           psurf=psurf+dtphys*dpsurf(1)   ! surface pressure change
977           DO ilevel=1,nlevel
978             plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
979           ENDDO
980           DO ilayer=1,nlayer
981             play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
982           ENDDO
983
984!       increment tracers
985        DO iq = 1, nq
986          DO ilayer=1,nlayer
987             q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
988          ENDDO
989        ENDDO
990      ENDDO   ! of idt=1,ndt ! end of time stepping loop
991
992      ! update the profiles files at the end of the run
993      write_prof=.false.
994      call getin("write_prof",write_prof)
995      IF (write_prof) then
996              DO iq = 1,nq
997                call writeprofile(nlayer,q(:,iq),noms(iq),qsurf)
998              ENDDO
999      ENDIF
1000
1001      write(*,*) "testphys1d: Everything is cool."
1002
1003      END
1004 
1005c***********************************************************************
1006c***********************************************************************
1007c     Dummy subroutines used only in 3D, but required to
1008c     compile testphys1d (to cleanly use writediagfi)
1009
1010      subroutine gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
1011
1012      IMPLICIT NONE
1013
1014      INTEGER im,jm,ngrid,nfield
1015      REAL pdyn(im,jm,nfield)
1016      REAL pfi(ngrid,nfield)
1017     
1018      if (ngrid.ne.1) then
1019        write(*,*) "gr_fi_dyn error: in 1D ngrid should be 1!!!"
1020        stop
1021      endif
1022     
1023      pdyn(1,1,1:nfield)=pfi(1,1:nfield)
1024     
1025      end
1026 
1027c***********************************************************************
1028c***********************************************************************
1029
Note: See TracBrowser for help on using the repository browser.