source: trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F @ 1529

Last change on this file since 1529 was 1529, checked in by emillour, 9 years ago

Generic GCM: Towards a cleaner separation between physics and dynamics

  • Got rid of references to "dimensions.h" from physics packages: use nbp_lon (=iim), nbp_lat (=jjp1) and nbp_lev from module mod_grid_phy_lmdz (in phy_common) instead.
  • Removed module "comhdiff_mod.F90", as it is only used by module surf_heat_transp_mod.F90, moved module variables there.
  • Addedin "surf_heat_transp_mod" local versions of some arrays and routines (from dyn3d) required to compute gradient, divergence, etc. on the global dynamics grid. As before, the slab ocean only works in serial.

EM
  • Property svn:executable set to *
File size: 29.0 KB
Line 
1      program rcm1d
2
3! to use  'getin'
4      use ioipsl_getincom, only: getin
5      use infotrac, only: nqtot, tname
6      use surfdat_h, only: albedodat, phisfi, dryness, watercaptag,
7     &                     zmea, zstd, zsig, zgam, zthe,
8     &                     emissiv, emisice, iceradius,
9     &                     dtemisice
10      use comdiurn_h, only: sinlat, coslat, sinlon, coslon
11!      use comsaison_h
12      use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat, volcapa
13      USE comgeomfi_h, only: lati, long, area
14      use phyredem, only: physdem0,physdem1
15      use comgeomphy, only: initcomgeomphy
16      use slab_ice_h, only: noceanmx
17      use planete_mod, only: apoastr,periastr,year_day,peri_day,
18     &         obliquit,nres,z0,lmixmin,emin_turb,coefvis,coefir,
19     &         timeperi,e_elips,p_elips
20      use comcstfi_mod, only: pi, cpp, rad, g, r,
21     &                        mugaz, rcp, omeg
22      use time_phylmdz_mod, only: daysec, dtphys, day_step, ecritphy,
23     &                            nday
24      use callkeys_mod, only: tracer,check_cpp_match,rings_shadow,
25     &                        specOLR,water,pceil,ok_slab_ocean
26      USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff
27      USE logic_mod, ONLY: hybrid,autozlevs
28      use inifis_mod, only: inifis
29      implicit none
30
31!==================================================================
32!     
33!     Purpose
34!     -------
35!     Run the physics package of the universal model in a 1D column.
36!     
37!     It can be compiled with a command like (e.g. for 25 layers):
38!                                  "makegcm -p std -d 25 rcm1d"
39!     It requires the files "callphys.def", "z2sig.def",
40!     "traceur.def", and "run.def" with a line "INCLUDEDEF=callphys.def"
41!
42!     Authors
43!     -------
44!     Frederic Hourdin
45!     R. Fournier
46!     F. Forget
47!     F. Montmessin (water ice added)
48!     R. Wordsworth
49!     B. Charnay
50!     A. Spiga
51!     J. Leconte (2012)
52!
53!==================================================================
54
55#include "dimensions.h"
56#include "paramet.h"
57!include "dimphys.h"
58#include "netcdf.inc"
59#include "comgeom.h"
60
61c --------------------------------------------------------------
62c  Declarations
63c --------------------------------------------------------------
64c
65      INTEGER unitstart      ! unite d'ecriture de "startfi"
66      INTEGER nlayer,nlevel,nsoil,ndt
67      INTEGER ilayer,ilevel,isoil,idt,iq
68      LOGICAl firstcall,lastcall
69c
70      INTEGER day0          ! date initial (sol ; =0 a Ls=0)
71      REAL day           ! date durant le run
72      REAL time             ! time (0<time<1 ; time=0.5 a midi)
73      REAL play(llm)   ! Pressure at the middle of the layers (Pa)
74      REAL plev(llm+1) ! intermediate pressure levels (pa)
75      REAL psurf,tsurf(1)     
76      REAL u(llm),v(llm)  ! zonal, meridional wind
77      REAL gru,grv   ! prescribed "geostrophic" background wind
78      REAL temp(llm)   ! temperature at the middle of the layers
79      REAL,ALLOCATABLE :: q(:,:) ! tracer mixing ratio (e.g. kg/kg)
80      REAL,ALLOCATABLE :: qsurf(:)      ! tracer surface budget (e.g. kg.m-2)
81      REAL tsoil(nsoilmx)   ! subsurface soik temperature (K)
82!      REAL co2ice           ! co2ice layer (kg.m-2) !not used anymore
83      integer :: i_co2_ice=0  ! tracer index of co2 ice
84      integer :: i_h2o_ice=0  ! tracer index of h2o ice
85      integer :: i_h2o_vap=0  ! tracer index of h2o vapor
86      REAL emis(1)             ! surface layer
87      REAL q2(llm+1)   ! Turbulent Kinetic Energy
88      REAL zlay(llm)   ! altitude estimee dans les couches (km)
89
90c    Physical and dynamical tandencies (e.g.  m.s-2, K/s, Pa/s)
91      REAL du(llm),dv(llm),dtemp(llm)
92      REAL dudyn(llm),dvdyn(llm),dtempdyn(llm)
93      REAL dpsurf   
94      REAL,ALLOCATABLE :: dq(:,:)
95      REAL,ALLOCATABLE :: dqdyn(:,:)
96
97c   Various intermediate variables
98!      INTEGER thermo
99      REAL zls
100      REAL phi(llm),h(llm),s(llm)
101      REAL pks, ptif, w(llm)
102      INTEGER ierr, aslun
103      REAL tmp1(0:llm),tmp2(0:llm)
104      Logical  tracerdyn
105      integer :: nq !=1 ! number of tracers
106 
107      character*2 str2
108      character (len=7) :: str7
109
110      logical oldcompare, earthhack,saveprofile
111
112!     added by RW for zlay computation
113      real Hscale, Hmax, rho, dz
114
115!     added by RW for autozlevs computation
116      real nu, xx, pMIN, zlev, Htop
117      real logplevs(llm)
118
119!     added by BC
120      REAL cloudfrac(1,llm)
121      REAL hice(1),totcloudfrac(1)
122      REAL qzero1D   !initial water amount on the ground
123
124!     added by BC for ocean
125      real rnat(1)
126      REAL tslab(1,noceanmx),tsea_ice(1),sea_ice(1)
127      real pctsrf_sic(1)
128
129
130
131!     added by AS to avoid the use of adv trac common
132      character*20,allocatable :: nametrac(:)   ! name of the tracer (no need for adv trac common)
133
134      real :: latitude(1), longitude(1)
135
136c=======================================================================
137c INITIALISATION
138c=======================================================================
139! initialize "serial/parallel" related stuff
140      CALL init_phys_lmdz(iim,jjm+1,llm,1,(/(jjm-1)*iim+2/))
141      call initcomgeomphy
142
143      !! those are defined in surfdat_h.F90
144      IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(1))
145      IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(1))
146      IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(1))
147      IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(1))
148      IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(1))
149      IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(1))
150      IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(1))
151      IF (.not. ALLOCATED(dryness)) ALLOCATE(dryness(1))
152      IF (.not. ALLOCATED(watercaptag)) ALLOCATE(watercaptag(1))
153      !! those are defined in comdiurn_h.F90
154      IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(1))
155      IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(1))
156      IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(1))
157      IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(1))
158      !! those are defined in comsoil_h.F90
159      IF (.not.ALLOCATED(layer)) ALLOCATE(layer(nsoilmx))
160      IF (.not.ALLOCATED(mlayer)) ALLOCATE(mlayer(0:nsoilmx-1))
161      IF (.not.ALLOCATED(inertiedat)) ALLOCATE(inertiedat(1,nsoilmx))
162      !! those are defined in comgeomfi_h
163      IF (.not. ALLOCATED(lati)) ALLOCATE(lati(1))
164      IF (.not. ALLOCATED(long)) ALLOCATE(long(1))
165      IF (.not. ALLOCATED(area)) ALLOCATE(area(1))
166
167
168      saveprofile=.false.
169      saveprofile=.true.
170
171c ----------------------------------------
172c  Default values  (corresponding to Mars)
173c ----------------------------------------
174
175      pi=2.E+0*asin(1.E+0)
176
177c     Parametres Couche limite et Turbulence
178c     --------------------------------------
179      z0 =  1.e-2                ! surface roughness (m) ~0.01
180      emin_turb = 1.e-6          ! energie minimale ~1.e-8
181      lmixmin = 30               ! longueur de melange ~100
182 
183c     propriete optiques des calottes et emissivite du sol
184c     ----------------------------------------------------
185      emissiv= 0.95              ! Emissivite du sol martien ~.95
186      emisice(1)=0.95            ! Emissivite calotte nord
187      emisice(2)=0.95            ! Emissivite calotte sud 
188
189      iceradius(1) = 100.e-6     ! mean scat radius of CO2 snow (north)
190      iceradius(2) = 100.e-6     ! mean scat radius of CO2 snow (south)
191      dtemisice(1) = 2.          ! time scale for snow metamorphism (north)
192      dtemisice(2) = 2.          ! time scale for snow metamorphism (south
193      hybrid=.false.
194
195c ------------------------------------------------------
196c  Load parameters from "run.def" and "gases.def"
197c ------------------------------------------------------
198
199! check if 'rcm1d.def' file is around
200      open(90,file='rcm1d.def',status='old',form='formatted',
201     &     iostat=ierr)
202      if (ierr.ne.0) then
203        write(*,*) 'Cannot find required file "rcm1d.def"'
204        write(*,*) 'which should contain some input parameters'
205        write(*,*) ' ... might as well stop here ...'
206        stop
207      else
208        close(90)
209      endif
210
211! now, run.def is needed anyway. so we create a dummy temporary one
212! for ioipsl to work. if a run.def is already here, stop the
213! program and ask the user to do a bit of cleaning
214      open(90,file='run.def',status='old',form='formatted',
215     &     iostat=ierr)
216      if (ierr.eq.0) then
217        close(90)
218        write(*,*) 'There is already a run.def file.'
219        write(*,*) 'This is not compatible with 1D runs.'
220        write(*,*) 'Please remove the file and restart the run.'
221        write(*,*) 'Runtime parameters are supposed to be in rcm1d.def'
222        stop
223      else
224        call system('touch run.def')
225        call system("echo 'INCLUDEDEF=callphys.def' >> run.def")
226        call system("echo 'INCLUDEDEF=rcm1d.def' >> run.def")
227      endif
228
229! check if we are going to run with or without tracers
230      write(*,*) "Run with or without tracer transport ?"
231      tracer=.false. ! default value
232      call getin("tracer",tracer)
233      write(*,*) " tracer = ",tracer
234
235! OK. now that run.def has been read once -- any variable is in memory.
236! so we can dump the dummy run.def
237      call system("rm -rf run.def")
238
239! while we're at it, check if there is a 'traceur.def' file
240! and preocess it, if necessary. Otherwise initialize tracer names
241      if (tracer) then
242      ! load tracer names from file 'traceur.def'
243        open(90,file='traceur.def',status='old',form='formatted',
244     &       iostat=ierr)
245        if (ierr.eq.0) then
246          write(*,*) "rcm1d: Reading file traceur.def"
247          ! read number of tracers:
248          read(90,*,iostat=ierr) nq
249          nqtot=nq ! set value of nqtot (in infotrac module) as nq
250          if (ierr.ne.0) then
251            write(*,*) "rcm1d: error reading number of tracers"
252            write(*,*) "   (first line of traceur.def) "
253            stop
254          endif
255          if (nq>0) then
256            allocate(tname(nq))
257            allocate(q(llm,nq))
258            allocate(qsurf(nq))
259            allocate(dq(llm,nq))
260            allocate(dqdyn(llm,nq))
261          else
262            write(*,*) "rcm1d: Error. You set tracer=.true."
263            write(*,*) "       but # of tracers in traceur.def is ",nq
264            stop
265          endif
266       
267          do iq=1,nq
268          ! minimal version, just read in the tracer names, 1 per line
269            read(90,*,iostat=ierr) tname(iq)
270            if (ierr.ne.0) then
271              write(*,*) 'rcm1d: error reading tracer names...'
272              stop
273            endif
274          enddo !of do iq=1,nq
275! check for co2_ice / h2o_ice tracers:
276         i_co2_ice=0
277         i_h2o_ice=0
278         i_h2o_vap=0
279         do iq=1,nq
280           if (tname(iq)=="co2_ice") then
281             i_co2_ice=iq
282           elseif (tname(iq)=="h2o_ice") then
283             i_h2o_ice=iq
284           elseif (tname(iq)=="h2o_vap") then
285             i_h2o_vap=iq
286           endif
287         enddo
288        else
289          write(*,*) 'Cannot find required file "traceur.def"'
290          write(*,*) ' If you want to run with tracers, I need it'
291          write(*,*) ' ... might as well stop here ...'
292          stop
293        endif
294        close(90)
295
296      else ! of if (tracer)
297        nqtot=0
298        nq=0
299        ! still, make allocations for 1 dummy tracer
300        allocate(tname(1))
301        allocate(qsurf(1))
302        allocate(q(llm,1))
303        allocate(dq(llm,1))
304     
305       ! Check that tracer boolean is compliant with number of tracers
306       ! -- otherwise there is an error (and more generally we have to be consistent)
307       if (nq .ge. 1) then
308          write(*,*) "------------------------------"
309          write(*,*) "rcm1d: You set tracer=.false."
310          write(*,*) " But set number of tracers to ",nq
311          write(*,*) " > If you want tracers, set tracer=.true."
312          write(*,*) "------------------------------"
313          stop
314       endif
315      endif ! of if (tracer)
316
317!!! We have to check that check_cpp_match is F for 1D computations
318!!! We think this check is better than make a particular case for 1D in inifis or calc_cpp_mugaz
319      check_cpp_match = .false.
320      call getin("check_cpp_match",check_cpp_match)
321      if (check_cpp_match) then
322          print*,"In 1D modeling, check_cpp_match is supposed to be F"
323          print*,"Please correct callphys.def"
324          stop
325      endif
326
327!!! GEOGRAPHICAL INITIALIZATIONS
328     !!! AREA. useless in 1D
329      area(1)=1.E+0
330      aire(1)=area(1) !JL+EM to have access to the area in the diagfi.nc files. area in comgeomfi.h and aire in comgeom.h
331     !!! GEOPOTENTIAL. useless in 1D because control by surface pressure
332      phisfi(1)=0.E+0
333     !!! LATITUDE. can be set.
334      latitude=0 ! default value for latitude
335      PRINT *,'latitude (in degrees) ?'
336      call getin("latitude",latitude)
337      write(*,*) " latitude = ",latitude
338      latitude=latitude*pi/180.E+0
339     !!! LONGITUDE. useless in 1D.
340      longitude=0.E+0
341      longitude=longitude*pi/180.E+0
342
343
344!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
345!!!! PLANETARY CONSTANTS !!!!
346!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
347
348      g = -99999.
349      PRINT *,'GRAVITY in m s-2 ?'
350      call getin("g",g)
351      IF (g.eq.-99999.) THEN
352          PRINT *,"STOP. I NEED g IN RCM1D.DEF."
353          STOP
354      ELSE
355          PRINT *,"--> g = ",g
356      ENDIF
357
358      rad = -99999.
359      PRINT *,'PLANETARY RADIUS in m ?'
360      call getin("rad",rad)
361      ! Planetary  radius is needed to compute shadow of the rings
362      IF (rad.eq.-99999. .and. rings_shadow .eqv. .true.) THEN
363          PRINT *,"STOP. I NEED rad IN RCM1D.DEF."
364          STOP
365      ELSE
366          PRINT *,"--> rad = ",rad
367      ENDIF
368
369      daysec = -99999.
370      PRINT *,'LENGTH OF A DAY in s ?'
371      call getin("daysec",daysec)
372      IF (daysec.eq.-99999.) THEN
373          PRINT *,"STOP. I NEED daysec IN RCM1D.DEF."
374          STOP
375      ELSE
376          PRINT *,"--> daysec = ",daysec
377      ENDIF
378      omeg=4.*asin(1.)/(daysec)
379      PRINT *,"OK. FROM THIS I WORKED OUT:"
380      PRINT *,"--> omeg = ",omeg
381
382      year_day = -99999.
383      PRINT *,'LENGTH OF A YEAR in days ?'
384      call getin("year_day",year_day)
385      IF (year_day.eq.-99999.) THEN
386          PRINT *,"STOP. I NEED year_day IN RCM1D.DEF."
387          STOP
388      ELSE
389          PRINT *,"--> year_day = ",year_day
390      ENDIF
391
392      periastr = -99999.
393      PRINT *,'MIN DIST STAR-PLANET in AU ?'
394      call getin("periastr",periastr)
395      IF (periastr.eq.-99999.) THEN
396          PRINT *,"STOP. I NEED periastr IN RCM1D.DEF."
397          STOP
398      ELSE
399          PRINT *,"--> periastr = ",periastr
400      ENDIF
401
402      apoastr = -99999.
403      PRINT *,'MAX DIST STAR-PLANET in AU ?'
404      call getin("apoastr",apoastr)
405      IF (apoastr.eq.-99999.) THEN
406          PRINT *,"STOP. I NEED apoastr IN RCM1D.DEF."
407          STOP
408      ELSE
409          PRINT *,"--> apoastr = ",apoastr
410      ENDIF
411
412      peri_day = -99999.
413      PRINT *,'DATE OF PERIASTRON in days ?'
414      call getin("peri_day",peri_day)
415      IF (peri_day.eq.-99999.) THEN
416          PRINT *,"STOP. I NEED peri_day IN RCM1D.DEF."
417          STOP
418      ELSE IF (peri_day.gt.year_day) THEN
419          PRINT *,"STOP. peri_day.gt.year_day"
420          STOP
421      ELSE 
422          PRINT *,"--> peri_day = ", peri_day
423      ENDIF
424
425      obliquit = -99999.
426      PRINT *,'OBLIQUITY in deg ?'
427      call getin("obliquit",obliquit)
428      IF (obliquit.eq.-99999.) THEN
429          PRINT *,"STOP. I NEED obliquit IN RCM1D.DEF."
430          STOP
431      ELSE
432          PRINT *,"--> obliquit = ",obliquit
433      ENDIF
434
435      psurf = -99999.
436      PRINT *,'SURFACE PRESSURE in Pa ?'
437      call getin("psurf",psurf)
438      IF (psurf.eq.-99999.) THEN
439          PRINT *,"STOP. I NEED psurf IN RCM1D.DEF."
440          STOP
441      ELSE
442          PRINT *,"--> psurf = ",psurf
443      ENDIF
444      !! we need reference pressures
445      pa=psurf/30.
446      preff=psurf ! these values are not needed in 1D  (are you sure JL12)
447
448!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
449!!!! END PLANETARY CONSTANTS !!!!
450!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
451
452c  Date et heure locale du debut du run
453c  ------------------------------------
454c    Date (en sols depuis le solstice de printemps) du debut du run
455      day0 = 0 ! default value for day0
456      write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?'
457      call getin("day0",day0)
458      day=float(day0)
459      write(*,*) " day0 = ",day0
460c  Heure de demarrage
461      time=0 ! default value for time
462      write(*,*)'Initial local time (in hours, between 0 and 24)?'
463      call getin("time",time)
464      write(*,*)" time = ",time
465      time=time/24.E+0 ! convert time (hours) to fraction of sol
466
467
468c  Discretisation (Definition de la grille et des pas de temps)
469c  --------------
470c
471      nlayer=llm
472      nlevel=nlayer+1
473      nsoil=nsoilmx
474
475      day_step=48 ! default value for day_step
476      PRINT *,'Number of time steps per sol ?'
477      call getin("day_step",day_step)
478      write(*,*) " day_step = ",day_step
479
480       
481      ecritphy=day_step ! default value for ecritphy
482      PRINT *,'Nunber of steps between writediagfi ?'
483      call getin("ecritphy",ecritphy)
484      write(*,*) " ecritphy = ",ecritphy
485
486      ndt=10 ! default value for ndt
487      PRINT *,'Number of sols to run ?'
488      call getin("ndt",ndt)
489      write(*,*) " ndt = ",ndt
490      nday=ndt
491
492      ndt=ndt*day_step     
493      dtphys=daysec/day_step 
494      write(*,*)"-------------------------------------"
495      write(*,*)"-------------------------------------"
496      write(*,*)"--> Physical timestep is ",dtphys
497      write(*,*)"-------------------------------------"
498      write(*,*)"-------------------------------------"
499
500!!! CALL INIFIS
501!!! - read callphys.def
502!!! - calculate sine and cosine of longitude and latitude
503!!! - calculate mugaz and cp
504!!! NB: some operations are being done dummily in inifis in 1D
505!!! - physical constants: nevermind, things are done allright below
506!!! - physical frequency: nevermind, in inifis this is a simple print
507      cpp=1.d-7 !JL because we divide by cpp in inifis, there may be a more elegant solution
508      CALL inifis(1,llm,nq,day0,daysec,nday,dtphys,
509     .            latitude,longitude,area,rad,g,r,cpp)
510!!! We check everything is OK.
511      PRINT *,"CHECK"
512      PRINT *,"--> mugaz = ",mugaz
513      PRINT *,"--> cpp = ",cpp
514      r = 8.314511E+0 * 1000.E+0 / mugaz
515      rcp = r / cpp
516
517c output spectrum?
518      write(*,*)"Output spectral OLR?"
519      specOLR=.false.
520      call getin("specOLR",specOLR)
521      write(*,*)" specOLR = ",specOLR
522
523c
524c  pour le schema d'ondes de gravite
525c  ---------------------------------
526      zmea(1)=0.E+0
527      zstd(1)=0.E+0
528      zsig(1)=0.E+0
529      zgam(1)=0.E+0
530      zthe(1)=0.E+0
531
532c    Initialisation des traceurs
533c    ---------------------------
534
535      DO iq=1,nq
536        DO ilayer=1,nlayer
537           q(ilayer,iq) = 0.
538        ENDDO
539      ENDDO
540     
541      DO iq=1,nq
542        qsurf(iq) = 0.
543      ENDDO
544
545      if (water) then
546         qzero1D               = 0.0
547         qsurf(i_h2o_vap) = qzero1D
548      endif
549
550c   Initialisation pour prendre en compte les vents en 1-D
551c   ------------------------------------------------------
552      ptif=2.E+0*omeg*sinlat(1)
553 
554
555c    vent geostrophique
556      gru=10. ! default value for gru
557      PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
558      call getin("u",gru)
559      write(*,*) " u = ",gru
560      grv=0. !default value for grv
561      PRINT *,'meridional northward component of the geostrophic',
562     &' wind (m/s) ?'
563      call getin("v",grv)
564      write(*,*) " v = ",grv
565
566c     Initialisation des vents  au premier pas de temps
567      DO ilayer=1,nlayer
568         u(ilayer)=gru
569         v(ilayer)=grv
570      ENDDO
571
572c     energie cinetique turbulente
573      DO ilevel=1,nlevel
574         q2(ilevel)=0.E+0
575      ENDDO
576
577c  emissivity / surface co2 ice ( + h2o ice??)
578c  -------------------------------------------
579      emis(1)=emissiv ! default value for emissivity
580      PRINT *,'Emissivity of bare ground ?'
581      call getin("emis",emis(1))
582      write(*,*) " emis = ",emis(1)
583      emissiv=emis(1) ! we do this so that condense_co2 sets things to the right
584                   ! value if there is no snow
585
586      if(i_co2_ice.gt.0)then
587         qsurf(i_co2_ice)=0 ! default value for co2ice
588         print*,'Initial CO2 ice on the surface (kg.m-2)'
589         call getin("co2ice",qsurf(i_co2_ice))
590         write(*,*) " co2ice = ",qsurf(i_co2_ice)
591         IF (qsurf(i_co2_ice).ge.1.E+0) THEN
592            ! if we have some CO2 ice on the surface, change emissivity
593            if (lati(1).ge.0) then ! northern hemisphere
594              emis(1)=emisice(1)
595            else ! southern hemisphere
596              emis(1)=emisice(2)
597            endif
598         ENDIF
599      endif
600
601c  calcul des pressions et altitudes en utilisant les niveaux sigma
602c  ----------------------------------------------------------------
603
604c    Vertical Coordinates
605c    """"""""""""""""""""
606      hybrid=.true.
607      PRINT *,'Hybrid coordinates ?'
608      call getin("hybrid",hybrid)
609      write(*,*) " hybrid = ", hybrid
610
611
612      autozlevs=.false.
613      PRINT *,'Auto-discretise vertical levels ?'
614      call getin("autozlevs",autozlevs)
615      write(*,*) " autozlevs = ", autozlevs
616
617      pceil = psurf / 1000.0 ! Pascals
618      PRINT *,'Ceiling pressure (Pa) ?'
619      call getin("pceil",pceil)
620      write(*,*) " pceil = ", pceil
621
622! Test of incompatibility:
623! if autozlevs used, cannot have hybrid too
624
625      if (autozlevs.and.hybrid) then
626         print*,'Cannot use autozlevs and hybrid together!'
627         call abort
628      endif
629
630      if(autozlevs)then
631           
632         open(91,file="z2sig.def",form='formatted')
633         read(91,*) Hscale
634         DO ilayer=1,nlayer-2
635            read(91,*) Hmax
636         enddo
637         close(91)
638 
639           
640         print*,'Hmax = ',Hmax,' km'
641         print*,'Auto-shifting Hscale to:'
642!     Hscale = Hmax / log(psurf/100.0)
643         Hscale = Hmax / log(psurf/pceil)
644         print*,'Hscale = ',Hscale,' km'
645         
646! none of this matters if we dont care about zlay
647         
648      endif
649
650      call disvert
651
652         if(.not.autozlevs)then
653            ! we want only the scale height from z2sig, in order to compute zlay
654            open(91,file="z2sig.def",form='formatted')
655            read(91,*) Hscale
656            close(91)
657         endif
658
659!         if(autozlevs)then
660!            open(94,file="Hscale.temp",form='formatted')
661!            read(94,*) Hscale
662!            close(94)
663!         endif
664
665         DO ilevel=1,nlevel
666            plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
667         ENDDO
668         
669         DO ilayer=1,nlayer
670            play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
671         ENDDO
672         
673
674
675         DO ilayer=1,nlayer
676!     zlay(ilayer)=-300.E+0 *r*log(play(ilayer)/plev(1))
677!     &   /g
678            zlay(ilayer)=-1000.0*Hscale*log(play(ilayer)/plev(1))
679         ENDDO
680
681!      endif
682
683c  profil de temperature au premier appel
684c  --------------------------------------
685      pks=psurf**rcp
686
687c altitude en km dans profile: on divise zlay par 1000
688      tmp1(0)=0.E+0
689      DO ilayer=1,nlayer
690        tmp1(ilayer)=zlay(ilayer)/1000.E+0
691      ENDDO
692      call profile(nlayer+1,tmp1,tmp2)
693
694      tsurf(1)=tmp2(0)
695      DO ilayer=1,nlayer
696        temp(ilayer)=tmp2(ilayer)
697      ENDDO
698      print*,"check"
699      PRINT*,"INPUT SURFACE TEMPERATURE",tsurf(1)
700      PRINT*,"INPUT TEMPERATURE PROFILE",temp
701
702c  Initialisation albedo / inertie du sol
703c  --------------------------------------
704      albedodat(1)=0.2 ! default value for albedodat
705      PRINT *,'Albedo of bare ground ?'
706      call getin("albedo",albedodat(1))
707      write(*,*) " albedo = ",albedodat(1)
708
709      inertiedat(1,1)=400 ! default value for inertiedat
710      PRINT *,'Soil thermal inertia (SI) ?'
711      call getin("inertia",inertiedat(1,1))
712      write(*,*) " inertia = ",inertiedat(1,1)
713
714! Initialize soil properties and temperature
715! ------------------------------------------
716      volcapa=1.e6 ! volumetric heat capacity
717      DO isoil=1,nsoil
718         inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
719         tsoil(isoil)=tsurf(1)  ! soil temperature
720      ENDDO
721
722! Initialize depths
723! -----------------
724      do isoil=0,nsoil-1
725        mlayer(isoil)=2.e-4*(2.**(isoil-0.5)) ! mid-layer depth
726      enddo
727      do isoil=1,nsoil
728        layer(isoil)=2.e-4*(2.**(isoil-1)) ! layer depth
729      enddo
730
731
732! Initialize slab ocean
733! -----------------
734      rnat=1. ! default value for rnat
735      if(inertiedat(1,1).GE.10000.)then
736         rnat=0.
737      endif
738      if(ok_slab_ocean)then
739      rnat=0.
740      tslab(1,1)=tsurf(1)
741      tslab(1,2)=tsurf(1)
742      tsea_ice=tsurf
743      pctsrf_sic=0.
744      sea_ice=0.
745      endif
746
747
748
749c  Write a "startfi" file
750c  --------------------
751c  This file will be read during the first call to "physiq".
752c  It is needed to transfert physics variables to "physiq"...
753
754      call physdem0("startfi.nc",long,lati,nsoilmx,1,llm,nq,
755     &              dtphys,real(day0),time,area,
756     &              albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
757      call physdem1("startfi.nc",nsoilmx,1,llm,nq,
758     &                dtphys,time,
759     &                tsurf,tsoil,emis,q2,qsurf,
760     &                cloudfrac,totcloudfrac,hice,
761     &                rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
762
763!      call physdem1(1,"startfi.nc",long,lati,nsoilmx,nq,
764!     &              dtphys,float(day0),
765!     &              time,tsurf,tsoil,emis,q2,qsurf,
766!     &              area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
767!     &              cloudfrac,totcloudfrac,hice,nametrac)
768
769c=======================================================================
770c  BOUCLE TEMPORELLE DU MODELE 1D
771c=======================================================================
772
773      firstcall=.true.
774      lastcall=.false.
775
776      DO idt=1,ndt
777        IF (idt.eq.ndt) then       !test
778         lastcall=.true.
779         call stellarlong(day*1.0,zls)
780!         write(103,*) 'Ls=',zls*180./pi
781!         write(103,*) 'Lat=', lati(1)*180./pi
782!         write(103,*) 'RunEnd - Atmos. Temp. File'
783!         write(103,*) 'RunEnd - Atmos. Temp. File'
784!         write(104,*) 'Ls=',zls*180./pi
785!         write(104,*) 'Lat=', lati(1)
786!         write(104,*) 'RunEnd - Atmos. Temp. File'
787        ENDIF
788
789c    calcul du geopotentiel
790c     ~~~~~~~~~~~~~~~~~~~~~
791
792
793      DO ilayer=1,nlayer
794
795!              if(autozlevs)then
796!                 s(ilayer)=(play(ilayer)/psurf)**rcp
797!              else
798          s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
799!              endif
800              !s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
801          h(ilayer)=cpp*temp(ilayer)/(pks*s(ilayer))
802       ENDDO
803
804!      DO ilayer=1,nlayer
805!        s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
806!        h(ilayer)=cpp*temp(ilayer)/(pks*s(ilayer))
807!      ENDDO
808      phi(1)=pks*h(1)*(1.E+0-s(1))
809      DO ilayer=2,nlayer
810         phi(ilayer)=phi(ilayer-1)+
811     &               pks*(h(ilayer-1)+h(ilayer))*.5E+0
812     &                  *(s(ilayer-1)-s(ilayer))
813
814      ENDDO
815
816c       appel de la physique
817c       --------------------
818
819
820      CALL physiq (1,llm,nq,
821     .     tname,
822     ,     firstcall,lastcall,
823     ,     day,time,dtphys,
824     ,     plev,play,phi,
825     ,     u, v,temp, q, 
826     ,     w,
827C - sorties
828     s     du, dv, dtemp, dq,dpsurf,tracerdyn)
829
830
831c       evolution du vent : modele 1D
832c       -----------------------------
833 
834c       la physique calcule les derivees temporelles de u et v.
835c       on y rajoute betement un effet Coriolis.
836c
837c       DO ilayer=1,nlayer
838c          du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
839c          dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
840c       ENDDO
841
842c       Pour certain test : pas de coriolis a l'equateur
843c       if(lati(1).eq.0.) then
844          DO ilayer=1,nlayer
845             du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
846             dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
847          ENDDO
848c       end if
849c     
850c
851c       Calcul du temps au pas de temps suivant
852c       ---------------------------------------
853        firstcall=.false.
854        time=time+dtphys/daysec
855        IF (time.gt.1.E+0) then
856            time=time-1.E+0
857            day=day+1
858        ENDIF
859
860c       calcul des vitesses et temperature au pas de temps suivant
861c       ----------------------------------------------------------
862
863        DO ilayer=1,nlayer
864           u(ilayer)=u(ilayer)+dtphys*du(ilayer)
865           v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
866           temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
867        ENDDO
868
869c       calcul des pressions au pas de temps suivant
870c       ----------------------------------------------------------
871
872           psurf=psurf+dtphys*dpsurf   ! evolution de la pression de surface
873           DO ilevel=1,nlevel
874              plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
875           ENDDO
876           DO ilayer=1,nlayer
877                 play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
878           ENDDO
879
880c       calcul traceur au pas de temps suivant
881c       --------------------------------------
882
883        DO iq = 1, nq
884          DO ilayer=1,nlayer
885             q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
886          ENDDO
887        END DO
888
889c    ========================================================
890c    GESTION DES SORTIE
891c    ========================================================
892      if(saveprofile)then
893         OPEN(12,file='profile.out',form='formatted')
894         write(12,*) tsurf
895         DO ilayer=1,llm
896            write(12,*) temp(ilayer) !, play(ilayer) !AS12 only temp so that iprofile=8 can be used
897         ENDDO
898         CLOSE(12)
899      endif
900
901
902      ENDDO   ! fin de la boucle temporelle
903
904      write(*,*) "rcm1d: Everything is cool."
905
906c    ========================================================
907      end                       !rcm1d
908 
909
Note: See TracBrowser for help on using the repository browser.