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

Last change on this file since 1271 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

  • Remove all reference/use of nlayermx and dimphys.h
  • Made use of automatic arrays whenever arrays are needed with dimension nlayer
  • Remove lots of obsolete reference to dimensions.h
  • Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

  • Deleted param.h and put contents into module param_v4_h
  • Adapted testphys1d, newstart, etc...
  • Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

  • Optimized computations in paramfoto_compact (twice less dlog10 calculations)
  • Checked consistency before/after modification in debug mode
  • Checked performance is not impacted (same as before)
File size: 27.4 KB
Line 
1
2      PROGRAM testphys1d
3! to use  'getin'
4      USE ioipsl_getincom, only: getin
5      use infotrac, only: nqtot, tname
6      use comsoil_h, only: volcapa, layer, mlayer, inertiedat, nsoilmx
7      use comgeomfi_h, only: lati, long, area, sinlat, ini_fillgeom
8      use surfdat_h, only: albedodat, z0_default, emissiv, emisice,
9     &                     albedice, iceradius, dtemisice, z0,
10     &                     zmea, zstd, zsig, zgam, zthe, phisfi,
11     &                     watercaptag
12      use slope_mod, only: theta_sl, psi_sl
13      use control_mod, only: day_step
14      use phyredem, only: physdem0,physdem1
15      use comgeomphy, only: initcomgeomphy
16      use planete_h, only: year_day, periheli, aphelie, peri_day,
17     &                     obliquit, emin_turb, lmixmin
18      use comcstfi_h, only: pi, rad, daysec, omeg, g, mugaz, rcp, r,
19     &                      cpp, dtphys
20      use dimradmars_mod, only: tauscaling,tauvis
21      IMPLICIT NONE
22
23c=======================================================================
24c   subject:
25c   --------
26c   PROGRAM useful to run physical part of the martian GCM in a 1D column
27c       
28c Can be compiled with a command like (e.g. for 25 layers)
29c  "makegcm -p mars -d 25 testphys1d"
30c It requires the files "testphys1d.def" "callphys.def"
31c   and a 'run.def' file (containing a "INCLUDEDEF=callphys.def" line)
32c      and a file describing the sigma layers (e.g. "z2sig.def")
33c
34c   author: Frederic Hourdin, R.Fournier,F.Forget
35c   -------
36c   
37c   update: 12/06/2003 including chemistry (S. Lebonnois)
38c                            and water ice (F. Montmessin)
39c
40c=======================================================================
41
42#include "dimensions.h"
43      integer, parameter :: ngrid = 1 !(2+(jjm-1)*iim - 1/jjm)
44      integer, parameter :: nlayer = llm
45!#include "dimradmars.h"
46!#include "comgeomfi.h"
47!#include "surfdat.h"
48!#include "slope.h"
49!#include "comsoil.h"
50!#include "comdiurn.h"
51#include "callkeys.h"
52!#include "comsaison.h"
53!#include "control.h"
54#include "comvert.h"
55#include "netcdf.inc"
56#include "comg1d.h"
57#include "logic.h"
58!#include "advtrac.h"
59
60c --------------------------------------------------------------
61c  Declarations
62c --------------------------------------------------------------
63c
64      INTEGER unitstart      ! unite d'ecriture de "startfi"
65      INTEGER nlayer,nlevel,nsoil,ndt
66      INTEGER ilayer,ilevel,isoil,idt,iq
67      LOGICAl firstcall,lastcall
68c
69      real,parameter :: odpref=610. ! DOD reference pressure (Pa)
70c
71      INTEGER day0          ! date initial (sol ; =0 a Ls=0)
72      REAL day           ! date durant le run
73      REAL time             ! time (0<time<1 ; time=0.5 a midi)
74      REAL play(nlayer)   ! Pressure at the middle of the layers (Pa)
75      REAL plev(nlayer+1) ! intermediate pressure levels (pa)
76      REAL psurf,tsurf(1)     
77      REAL u(nlayer),v(nlayer)  ! zonal, meridional wind
78      REAL gru,grv   ! prescribed "geostrophic" background wind
79      REAL temp(nlayer)   ! temperature at the middle of the layers
80      REAL,ALLOCATABLE :: q(:,:) ! tracer mixing ratio (e.g. kg/kg)
81      REAL,ALLOCATABLE :: qsurf(:) ! tracer surface budget (e.g. kg.m-2)
82      REAL tsoil(nsoilmx)   ! subsurface soik temperature (K)
83      REAL co2ice(1)        ! co2ice layer (kg.m-2)
84      REAL emis(1)          ! surface layer
85      REAL q2(nlayer+1)   ! Turbulent Kinetic Energy
86      REAL zlay(nlayer)   ! altitude estimee dans les couches (km)
87
88c    Physical and dynamical tandencies (e.g.  m.s-2, K/s, Pa/s)
89      REAL du(nlayer),dv(nlayer),dtemp(nlayer)
90      REAL dudyn(nlayer),dvdyn(nlayer),dtempdyn(nlayer)
91      REAL dpsurf   
92      REAL,ALLOCATABLE :: dq(:,:)
93      REAL,ALLOCATABLE :: dqdyn(:,:)
94
95c   Various intermediate variables
96      INTEGER thermo
97      REAL zls
98      REAL phi(nlayer),h(nlayer),s(nlayer)
99      REAL pks, ptif, w(nlayer)
100      REAL qtotinit,qtot
101      real,allocatable :: mqtot(:)
102      INTEGER ierr, aslun
103      REAL tmp1(0:nlayer),tmp2(0:nlayer)
104      Logical  tracerdyn
105      integer :: nq=1 ! number of tracers
106      real :: latitude, longitude
107
108      character*2 str2
109      character (len=7) :: str7
110      character(len=44) :: txt
111
112c=======================================================================
113
114c=======================================================================
115c INITIALISATION
116c=======================================================================
117! initialize "serial/parallel" related stuff
118!      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
119      CALL init_phys_lmdz(iim,jjm+1,llm,1,(/(jjm-1)*iim+2/))
120      call initcomgeomphy
121
122c ------------------------------------------------------
123c  Prescribed constants to be set here
124c ------------------------------------------------------
125
126      pi=2.E+0*asin(1.E+0)
127
128c     Mars planetary constants
129c     ----------------------------
130      rad=3397200.               ! mars radius (m)  ~3397200 m
131      daysec=88775.              ! length of a sol (s)  ~88775 s
132      omeg=4.*asin(1.)/(daysec)  ! rotation rate (rad.s-1)
133      g=3.72                     ! gravity (m.s-2) ~3.72 
134      mugaz=43.49                ! atmosphere mola mass (g.mol-1) ~43.49
135      rcp=.256793                ! = r/cp  ~0.256793
136      r= 8.314511E+0 *1000.E+0/mugaz
137      cpp= r/rcp
138      year_day = 669             ! lenght of year (sols) ~668.6
139      periheli = 206.66          ! minimum sun-mars distance (Mkm) ~206.66
140      aphelie = 249.22           ! maximum sun-mars distance (Mkm) ~249.22
141      peri_day =  485.           ! perihelion date (sols since N. Spring)
142      obliquit = 25.2            ! Obliquity (deg) ~25.2         
143 
144c     Planetary Boundary Layer and Turbulence parameters
145c     --------------------------------------------------
146      z0_default =  1.e-2        ! surface roughness (m) ~0.01
147      emin_turb = 1.e-6          ! minimal turbulent energy ~1.e-8
148      lmixmin = 30               ! mixing length ~100
149 
150c     cap properties and surface emissivities
151c     ----------------------------------------------------
152      emissiv= 0.95              ! Bare ground emissivity ~.95
153      emisice(1)=0.95            ! Northern cap emissivity
154      emisice(2)=0.95            ! Southern cap emisssivity
155      albedice(1)=0.5            ! Northern cap albedo
156      albedice(2)=0.5            ! Southern cap albedo
157      iceradius(1) = 100.e-6     ! mean scat radius of CO2 snow (north)
158      iceradius(2) = 100.e-6     ! mean scat radius of CO2 snow (south)
159      dtemisice(1) = 2.          ! time scale for snow metamorphism (north)
160      dtemisice(2) = 2.          ! time scale for snow metamorphism (south
161
162
163c ------------------------------------------------------
164c  Loading run parameters from "run.def" file
165c ------------------------------------------------------
166
167
168! check if 'run.def' file is around (otherwise reading parameters
169! from callphys.def via getin() routine won't work.
170      open(99,file='run.def',status='old',form='formatted',
171     &     iostat=ierr)
172      if (ierr.ne.0) then
173        write(*,*) 'Cannot find required file "run.def"'
174        write(*,*) '  (which should contain some input parameters'
175        write(*,*) '   along with the following line:'
176        write(*,*) '   INCLUDEDEF=callphys.def'
177        write(*,*) '   )'
178        write(*,*) ' ... might as well stop here ...'
179        stop
180      else
181        close(99)
182      endif
183
184! check if we are going to run with or without tracers
185      write(*,*) "Run with or without tracer transport ?"
186      tracer=.false. ! default value
187      call getin("tracer",tracer)
188      write(*,*) " tracer = ",tracer
189
190! while we're at it, check if there is a 'traceur.def' file
191! and process it.
192      if (tracer) then
193      ! load tracer names from file 'traceur.def'
194        open(90,file='traceur.def',status='old',form='formatted',
195     &       iostat=ierr)
196        if (ierr.ne.0) then
197          write(*,*) 'Cannot find required file "traceur.def"'
198          write(*,*) ' If you want to run with tracers, I need it'
199          write(*,*) ' ... might as well stop here ...'
200          stop
201        else
202          write(*,*) "testphys1d: Reading file traceur.def"
203          ! read number of tracers:
204          read(90,*,iostat=ierr) nq
205          nqtot=nq ! set value of nqtot (in infotrac module) as nq
206          if (ierr.ne.0) then
207            write(*,*) "testphys1d: error reading number of tracers"
208            write(*,*) "   (first line of traceur.def) "
209            stop
210          endif
211        endif
212        ! allocate arrays:
213        allocate(tname(nq))
214        allocate(q(nlayer,nq))
215        allocate(qsurf(nq))
216        allocate(dq(nlayer,nq))
217        allocate(dqdyn(nlayer,nq))
218        allocate(mqtot(nq))
219       
220        ! read tracer names from file traceur.def
221        do iq=1,nq
222          read(90,*,iostat=ierr) tname(iq)
223          if (ierr.ne.0) then
224            write(*,*) 'testphys1d: error reading tracer names...'
225            stop
226          endif
227        enddo
228        close(90)
229
230        ! initialize tracers here:
231        write(*,*) "testphys1d: initializing tracers"
232        q(:,:)=0 ! default, set everything to zero
233        qsurf(:)=0
234        ! "smarter" initialization of some tracers
235        ! (get values from "profile_*" files, if these are available)
236        do iq=1,nq
237          txt=""
238          write(txt,"(a)") tname(iq)
239          write(*,*)"  tracer:",trim(txt)
240          ! CO2
241          if (txt.eq."co2") then
242            q(:,iq)=0.95   ! kg /kg of atmosphere
243            qsurf(iq)=0. ! kg/m2 (not used for CO2)
244            ! even better, look for a "profile_co2" input file
245            open(91,file='profile_co2',status='old',
246     &       form='formatted',iostat=ierr)
247            if (ierr.eq.0) then
248              read(91,*) qsurf(iq)
249              do ilayer=1,nlayer
250                read(91,*) q(ilayer,iq)
251              enddo
252            endif
253            close(91)
254          endif ! of if (txt.eq."co2")
255          ! Allow for an initial profile of argon
256          ! Can also be used to introduce a decaying tracer
257          ! in the 1D (TBD) to study thermals
258          if (txt.eq."ar") then
259            !look for a "profile_ar" input file
260            open(91,file='profile_ar',status='old',
261     &       form='formatted',iostat=ierr)
262            if (ierr.eq.0) then
263              read(91,*) qsurf(iq)
264              do ilayer=1,nlayer
265                read(91,*) q(ilayer,iq)
266              enddo
267            else
268              write(*,*) "No profile_ar file!"
269            endif
270            close(91)
271          endif ! of if (txt.eq."ar")
272
273          ! WATER VAPOUR
274          if (txt.eq."h2o_vap") then
275            !look for a "profile_h2o_vap" input file
276            open(91,file='profile_h2o_vap',status='old',
277     &       form='formatted',iostat=ierr)
278            if (ierr.eq.0) then
279              read(91,*) qsurf(iq)
280              do ilayer=1,nlayer
281                read(91,*) q(ilayer,iq)
282              enddo
283            else
284              write(*,*) "No profile_h2o_vap file!"
285            endif
286            close(91)
287          endif ! of if (txt.eq."h2o_ice")
288          ! WATER ICE
289          if (txt.eq."h2o_ice") then
290            !look for a "profile_h2o_vap" input file
291            open(91,file='profile_h2o_ice',status='old',
292     &       form='formatted',iostat=ierr)
293            if (ierr.eq.0) then
294              read(91,*) qsurf(iq)
295              do ilayer=1,nlayer
296                read(91,*) q(ilayer,iq)
297              enddo
298            else
299              write(*,*) "No profile_h2o_ice file!"
300            endif
301            close(91)
302          endif ! of if (txt.eq."h2o_ice")
303          ! DUST
304          !if (txt(1:4).eq."dust") then
305          !  q(:,iq)=0.4    ! kg/kg of atmosphere
306          !  qsurf(iq)=100 ! kg/m2
307          !endif
308          ! DUST MMR
309          if (txt.eq."dust_mass") then
310            !look for a "profile_dust_mass" input file
311            open(91,file='profile_dust_mass',status='old',
312     &       form='formatted',iostat=ierr)
313            if (ierr.eq.0) then
314              read(91,*) qsurf(iq)
315              do ilayer=1,nlayer
316                read(91,*) q(ilayer,iq)
317!                write(*,*) "l=",ilayer," q(ilayer,iq)=",q(ilayer,iq)
318              enddo
319            else
320              write(*,*) "No profile_dust_mass file!"
321            endif
322            close(91)
323          endif ! of if (txt.eq."dust_mass")
324          ! DUST NUMBER
325          if (txt.eq."dust_number") then
326            !look for a "profile_dust_number" input file
327            open(91,file='profile_dust_number',status='old',
328     &       form='formatted',iostat=ierr)
329            if (ierr.eq.0) then
330              read(91,*) qsurf(iq)
331              do ilayer=1,nlayer
332                read(91,*) q(ilayer,iq)
333              enddo
334            else
335              write(*,*) "No profile_dust_number file!"
336            endif
337            close(91)
338          endif ! of if (txt.eq."dust_number")
339          ! NB: some more initializations (chemistry) is done later
340          ! CCN MASS
341          if (txt.eq."ccn_mass") then
342            !look for a "profile_ccn_mass" input file
343            open(91,file='profile_ccn_mass',status='old',
344     &       form='formatted',iostat=ierr)
345            if (ierr.eq.0) then
346              read(91,*) qsurf(iq)
347              do ilayer=1,nlayer
348                read(91,*) q(ilayer,iq)
349              enddo
350            else
351              write(*,*) "No profile_ccn_mass file!"
352            endif
353            close(91)
354          endif ! of if (txt.eq."ccn_mass")
355          ! CCN NUMBER
356          if (txt.eq."ccn_number") then
357            !look for a "profile_ccn_number" input file
358            open(91,file='profile_ccn_number',status='old',
359     &       form='formatted',iostat=ierr)
360            if (ierr.eq.0) then
361              read(91,*) qsurf(iq)
362              do ilayer=1,nlayer
363                read(91,*) q(ilayer,iq)
364              enddo
365            else
366              write(*,*) "No profile_ccn_number file!"
367            endif
368            close(91)
369          endif ! of if (txt.eq."ccn_number")
370        enddo ! of do iq=1,nq
371
372      else
373      ! we still need to set (dummy) tracer number and names for physdem1
374        nq=1
375        ! allocate arrays:
376        allocate(tname(nq))
377        allocate(q(nlayer,nq))
378        allocate(qsurf(nq))
379        allocate(dq(nlayer,nq))
380        allocate(dqdyn(nlayer,nq))
381        allocate(mqtot(nq))
382        do iq=1,nq
383          write(str7,'(a1,i2.2)')'q',iq
384          tname(iq)=str7
385        enddo
386      ! and just to be clean, also initialize tracers to zero for physdem1
387        q(:,:)=0
388        qsurf(:)=0     
389      endif ! of if (tracer)
390     
391      !write(*,*) "testphys1d q", q(1,:)
392      !write(*,*) "testphys1d qsurf", qsurf
393
394c  Date and local time at beginning of run
395c  ---------------------------------------
396c    Date (in sols since spring solstice) at beginning of run
397      day0 = 0 ! default value for day0
398      write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?'
399      call getin("day0",day0)
400      day=float(day0)
401      write(*,*) " day0 = ",day0
402c  Local time at beginning of run
403      time=0 ! default value for time
404      write(*,*)'Initial local time (in hours, between 0 and 24)?'
405      call getin("time",time)
406      write(*,*)" time = ",time
407      time=time/24.E+0 ! convert time (hours) to fraction of sol
408
409c  Discretization (Definition of grid and time steps)
410c  --------------
411c
412      nlevel=nlayer+1
413      nsoil=nsoilmx
414
415      day_step=48 ! default value for day_step
416      PRINT *,'Number of time steps per sol ?'
417      call getin("day_step",day_step)
418      write(*,*) " day_step = ",day_step
419
420      ndt=10 ! default value for ndt
421      PRINT *,'Number of sols to run ?'
422      call getin("ndt",ndt)
423      write(*,*) " ndt = ",ndt
424
425      ndt=ndt*day_step     
426      dtphys=daysec/day_step 
427
428c Imposed surface pressure
429c ------------------------------------
430c
431      psurf=610. ! default value for psurf
432      PRINT *,'Surface pressure (Pa) ?'
433      call getin("psurf",psurf)
434      write(*,*) " psurf = ",psurf
435c Reference pressures
436      pa=20.   ! transition pressure (for hybrid coord.)
437      preff=610.      ! reference surface pressure
438 
439c Aerosol properties
440c --------------------------------
441      tauvis=0.2 ! default value for tauvis (dust opacity)
442      write(*,'("Reference dust opacity at ",f4.0," Pa ?")')odpref
443      call getin("tauvis",tauvis)
444      write(*,*) " tauvis = ",tauvis
445
446c Orbital parameters
447c ------------------
448      print *,'Min. distance Sun-Mars (Mkm)?'
449      call getin("periheli",periheli)
450      write(*,*) " periheli = ",periheli
451
452      print *,'Max. distance Sun-Mars (Mkm)?'
453      call getin("aphelie",aphelie)
454      write(*,*) " aphelie = ",aphelie
455
456      print *,'Day of perihelion?'
457      call getin("periday",peri_day)
458      write(*,*) " periday = ",peri_day
459
460      print *,'Obliquity?'
461      call getin("obliquit",obliquit)
462      write(*,*) " obliquit = ",obliquit
463 
464c  latitude/longitude
465c  ------------------
466      latitude=0 ! default value for latitude
467      PRINT *,'latitude (in degrees) ?'
468      call getin("latitude",latitude)
469      write(*,*) " latitude = ",latitude
470      latitude=latitude*pi/180.E+0
471      longitude=0.E+0
472      longitude=longitude*pi/180.E+0
473
474!  some initializations (some of which have already been
475!  done above!) and loads parameters set in callphys.def
476!  and allocates some arrays
477!Mars possible matter with dtphys in input and include!!!
478      call phys_state_var_init(1,llm,nq,
479     .          daysec,dtphys,rad,g,r,cpp)
480      call ini_fillgeom(1,latitude,longitude,1.0)
481      call conf_phys(1,llm,nq)
482
483
484c  Initialize albedo / soil thermal inertia
485c  ----------------------------------------
486c
487      albedodat(1)=0.2 ! default value for albedodat
488      PRINT *,'Albedo of bare ground ?'
489      call getin("albedo",albedodat(1))
490      write(*,*) " albedo = ",albedodat(1)
491
492      inertiedat(1,1)=400 ! default value for inertiedat
493      PRINT *,'Soil thermal inertia (SI) ?'
494      call getin("inertia",inertiedat(1,1))
495      write(*,*) " inertia = ",inertiedat(1,1)
496
497      z0(1)=z0_default ! default value for roughness
498      write(*,*) 'Surface roughness length z0 (m)?'
499      call getin("z0",z0(1))
500      write(*,*) " z0 = ",z0(1)
501
502! Initialize local slope parameters (only matters if "callslope"
503! is .true. in callphys.def)
504      ! slope inclination angle (deg) 0: horizontal, 90: vertical
505      theta_sl(1)=0.0 ! default: no inclination
506      call getin("slope_inclination",theta_sl(1))
507      ! slope orientation (deg)
508      ! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
509      psi_sl(1)=0.0 ! default value
510      call getin("slope_orientation",psi_sl(1))
511     
512c
513c  for the gravity wave scheme
514c  ---------------------------------
515c
516      zmea(1)=0.E+0
517      zstd(1)=0.E+0
518      zsig(1)=0.E+0
519      zgam(1)=0.E+0
520      zthe(1)=0.E+0
521
522
523c   Specific initializations for "physiq"
524c   -------------------------------------
525c   mesh surface (not a very usefull quantity in 1D)
526      area(1)=1.E+0
527
528c   surface geopotential is not used (or useful) since in 1D
529c   everything is controled by surface pressure
530      phisfi(1)=0.E+0
531
532c   Initialization to take into account prescribed winds
533c   ------------------------------------------------------
534      ptif=2.E+0*omeg*sinlat(1)
535 
536c    geostrophic wind
537      gru=10. ! default value for gru
538      PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
539      call getin("u",gru)
540      write(*,*) " u = ",gru
541      grv=0. !default value for grv
542      PRINT *,'meridional northward component of the geostrophic',
543     &' wind (m/s) ?'
544      call getin("v",grv)
545      write(*,*) " v = ",grv
546
547c     Initialize winds  for first time step
548      DO ilayer=1,nlayer
549         u(ilayer)=gru
550         v(ilayer)=grv
551      ENDDO
552
553c     Initialize turbulente kinetic energy
554      DO ilevel=1,nlevel
555         q2(ilevel)=0.E+0
556      ENDDO
557
558c  CO2 ice on the surface
559c  -------------------
560      co2ice(1)=0.E+0 ! default value for co2ice
561      PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
562      call getin("co2ice",co2ice)
563      write(*,*) " co2ice = ",co2ice
564
565c
566c  emissivity
567c  ----------
568      emis=emissiv
569      IF (co2ice(1).eq.1.E+0) THEN
570         emis=emisice(1) ! northern hemisphere
571         IF(lati(1).LT.0) emis=emisice(2) ! southern hemisphere
572      ENDIF
573
574 
575
576c  Compute pressures and altitudes of atmospheric levels
577c  ----------------------------------------------------------------
578
579c    Vertical Coordinates
580c    """"""""""""""""""""
581      hybrid=.true.
582      PRINT *,'Hybrid coordinates ?'
583      call getin("hybrid",hybrid)
584      write(*,*) " hybrid = ", hybrid
585
586      CALL  disvert
587
588      DO ilevel=1,nlevel
589        plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
590      ENDDO
591
592      DO ilayer=1,nlayer
593        play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
594      ENDDO
595
596      DO ilayer=1,nlayer
597        zlay(ilayer)=-200.E+0 *r*log(play(ilayer)/plev(1))
598     &   /g
599      ENDDO
600
601
602c  Initialize temperature profile
603c  --------------------------------------
604      pks=psurf**rcp
605
606c altitude in km in profile: divide zlay by 1000
607      tmp1(0)=0.E+0
608      DO ilayer=1,nlayer
609        tmp1(ilayer)=zlay(ilayer)/1000.E+0
610      ENDDO
611
612      call profile(nlayer+1,tmp1,tmp2)
613
614      tsurf=tmp2(0)
615      DO ilayer=1,nlayer
616        temp(ilayer)=tmp2(ilayer)
617      ENDDO
618     
619
620
621! Initialize soil properties and temperature
622! ------------------------------------------
623      volcapa=1.e6 ! volumetric heat capacity
624      DO isoil=1,nsoil
625         inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
626         tsoil(isoil)=tsurf(1)  ! soil temperature
627      ENDDO
628
629! Initialize depths
630! -----------------
631      do isoil=0,nsoil-1
632        mlayer(isoil)=2.e-4*(2.**(isoil-0.5)) ! mid-layer depth
633      enddo
634      do isoil=1,nsoil
635        layer(isoil)=2.e-4*(2.**(isoil-1)) ! layer depth
636      enddo
637
638c    Initialize traceurs
639c    ---------------------------
640
641      if (photochem.or.callthermos) then
642         write(*,*) 'Initializing chemical species'
643         ! thermo=0: initialize over all atmospheric layers
644         thermo=0
645         call inichim_newstart(q,psurf,sig,nq,lati,long,area,
646     $   thermo,qsurf)
647      endif
648
649c Check if the surface is a water ice reservoir
650c --------------------------------------------------
651      watercaptag(1)=.false. ! Default: no water ice reservoir
652      print *,'Water ice cap on ground ?'
653      call getin("watercaptag",watercaptag)
654      write(*,*) " watercaptag = ",watercaptag
655     
656
657c    Initialization for GRADS outputs in "g1d.dat" and "g1d.ctl"
658c    ----------------------------------------------------------------
659c    (output done in "writeg1d", typically called by "physiq.F")
660
661        g1d_nlayer=nlayer
662        g1d_nomfich='g1d.dat'
663        g1d_unitfich=40
664        g1d_nomctl='g1d.ctl'
665        g1d_unitctl=41
666        g1d_premier=.true.
667        g2d_premier=.true.
668
669c  Write a "startfi" file
670c  --------------------
671c  This file will be read during the first call to "physiq".
672c  It is needed to transfert physics variables to "physiq"...
673
674      call physdem0("startfi.nc",long,lati,nsoilmx,ngrid,llm,nq,
675     .              dtphys,float(day0),time,area,
676     .              albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
677      call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,
678     .              dtphys,time,
679     .              tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling)
680
681c=======================================================================
682c  1D MODEL TIME STEPPING LOOP
683c=======================================================================
684c
685      firstcall=.true.
686      lastcall=.false.
687
688      DO idt=1,ndt
689c        IF (idt.eq.ndt) lastcall=.true.
690        IF (idt.eq.ndt-day_step-1) then       !test
691         lastcall=.true.
692         call solarlong(day*1.0,zls)
693         write(103,*) 'Ls=',zls*180./pi
694         write(103,*) 'Lat=', lati(1)*180./pi
695         write(103,*) 'Tau=', tauvis/odpref*psurf
696         write(103,*) 'RunEnd - Atmos. Temp. File'
697         write(103,*) 'RunEnd - Atmos. Temp. File'
698         write(104,*) 'Ls=',zls*180./pi
699         write(104,*) 'Lat=', lati(1)
700         write(104,*) 'Tau=', tauvis/odpref*psurf
701         write(104,*) 'RunEnd - Atmos. Temp. File'
702        ENDIF
703
704c     compute geopotential
705c     ~~~~~~~~~~~~~~~~~~~~~
706      DO ilayer=1,nlayer
707        s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
708        h(ilayer)=cpp*temp(ilayer)/(pks*s(ilayer))
709      ENDDO
710      phi(1)=pks*h(1)*(1.E+0-s(1))
711      DO ilayer=2,nlayer
712         phi(ilayer)=phi(ilayer-1)+
713     &               pks*(h(ilayer-1)+h(ilayer))*.5E+0
714     &                  *(s(ilayer-1)-s(ilayer))
715
716      ENDDO
717
718c       call physics
719c       --------------------
720!      write(*,*) "testphys1d avant q", q(1,:)
721      CALL physiq (1,llm,nq,
722     ,     firstcall,lastcall,
723     ,     day,time,dtphys,
724     ,     plev,play,phi,
725     ,     u, v,temp, q, 
726     ,     w,
727C - outputs
728     s     du, dv, dtemp, dq,dpsurf,tracerdyn)
729!      write(*,*) "testphys1d apres q", q(1,:)
730
731
732c       wind increment : specific for 1D
733c       --------------------------------
734 
735c       The physics compute the tendencies on u and v,
736c       here we just add Coriolos effect
737c
738c       DO ilayer=1,nlayer
739c          du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
740c          dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
741c       ENDDO
742
743c       For some tests : No coriolis force at equator
744c       if(lati(1).eq.0.) then
745          DO ilayer=1,nlayer
746             du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
747             dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
748          ENDDO
749c       end if
750c     
751c
752c       Compute time for next time step
753c       ---------------------------------------
754        firstcall=.false.
755        time=time+dtphys/daysec
756        IF (time.gt.1.E+0) then
757            time=time-1.E+0
758            day=day+1
759        ENDIF
760
761c       compute winds and temperature for next time step
762c       ----------------------------------------------------------
763
764        DO ilayer=1,nlayer
765           u(ilayer)=u(ilayer)+dtphys*du(ilayer)
766           v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
767           temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
768        ENDDO
769
770c       compute pressure for next time step
771c       ----------------------------------------------------------
772
773           psurf=psurf+dtphys*dpsurf   ! surface pressure change
774           DO ilevel=1,nlevel
775             plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
776           ENDDO
777           DO ilayer=1,nlayer
778             play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
779           ENDDO
780
781!       increment tracers
782        DO iq = 1, nq
783          DO ilayer=1,nlayer
784             q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
785          ENDDO
786        ENDDO
787
788      ENDDO   ! of idt=1,ndt ! end of time stepping loop
789
790c    ========================================================
791c    OUTPUTS
792c    ========================================================
793
794c    finalize and close grads files "g1d.dat" and "g1d.ctl"
795
796c        CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
797        CALL endg1d(1,nlayer,zlay/1000.,ndt)
798
799c    ========================================================
800      END
801 
802c***********************************************************************
803c***********************************************************************
804c     Dummy subroutines used only in 3D, but required to
805c     compile testphys1d (to cleanly use writediagfi)
806
807      subroutine gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
808
809      IMPLICIT NONE
810
811      INTEGER im,jm,ngrid,nfield
812      REAL pdyn(im,jm,nfield)
813      REAL pfi(ngrid,nfield)
814     
815      if (ngrid.ne.1) then
816        write(*,*) "gr_fi_dyn error: in 1D ngrid should be 1!!!"
817        stop
818      endif
819     
820      pdyn(1,1,1:nfield)=pfi(1,1:nfield)
821     
822      end
823 
824c***********************************************************************
825c***********************************************************************
826
827#include "../dyn3d/disvert.F"
828#include "../dyn3d/abort_gcm.F"
Note: See TracBrowser for help on using the repository browser.