source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/phys/kcm1d.f90 @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago
File size: 16.9 KB
Line 
1program kcm1d
2
3  use infotrac, only: nqtot
4  use radinc_h,      only: NAERKIND
5  use radcommon_h,   only: L_NSPECTI, L_NSPECTV, sigma
6  use watercommon_h, only: mH2O
7  use ioipsl_getincom, only: getin
8  use comsaison_h, only: mu0, fract, dist_star
9!  use control_mod
10
11  implicit none
12
13  !==================================================================
14  !     
15  !     Purpose
16  !     -------
17  !     Run the universal model radiative transfer once in a 1D column.
18  !     Useful for climate evolution studies etc.
19  !     
20  !     It can be compiled with a command like (e.g. for 25 layers):
21  !                                  "makegcm -p std -d 25 kcm1d"
22  !     It requires the files "callphys.def", "gases.def"
23  !     "traceur.def", and "run.def" with a line "INCLUDEDEF=callphys.def"
24  !
25  !     Authors
26  !     -------
27  !     R. Wordsworth
28  !
29  !==================================================================
30
31!-----------------------------------------------------------------------
32!   INCLUDE 'dimensions.h'
33!
34!   dimensions.h contient les dimensions du modele
35!   ndm est tel que iim=2**ndm
36!-----------------------------------------------------------------------
37
38      INTEGER iim,jjm,llm,ndm
39
40      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
41
42!-----------------------------------------------------------------------
43!-----------------------------------------------------------------------
44!   INCLUDE 'dimphys.h'
45
46! ngridmx : number of horizontal grid points
47! note: the -1/jjm term will be 0; unless jj=1
48      integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)   
49! nlayermx : number of atmospheric layers
50      integer, parameter :: nlayermx = llm 
51! nsoilmx : number of subterranean layers ! nsoilmx is now in comsoil_h
52      !integer, parameter :: nsoilmx = 4 ! for a test
53      !integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case
54      !integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
55!-----------------------------------------------------------------------
56
57!
58! For Fortran 77/Fortran 90 compliance always use line continuation
59! symbols '&' in columns 73 and 6
60!
61! Group commons according to their type for minimal performance impact
62
63      COMMON/callkeys_l/callrad,corrk,calldifv,UseTurbDiff,calladj      &
64     &   , co2cond,callsoil                                             &
65     &   , season,diurnal,tlocked,rings_shadow,lwrite                   &
66     &   , callstats,calleofdump                                        &
67     &   , enertest                                                     &
68     &   , callgasvis,continuum,H2Ocont_simple,graybody                 &
69     &   , radfixed                                                     &
70     &   , meanOLR, specOLR                                             &
71     &   , kastprof                                                     &
72     &   , nosurf, oblate                                               &     
73     &   , newtonian, testradtimes                                      &
74     &   , check_cpp_match, force_cpp                                   &
75     &   , rayleigh                                                     &
76     &   , stelbbody                                                    &
77     &   , nearco2cond                                                  &
78     &   , tracer, mass_redistrib, varactive, varfixed                  &
79     &   , sedimentation,water,watercond,waterrain                      &
80     &   , aeroco2,aeroh2o,aeroh2so4,aeroback2lay                       &
81     &   , aerofixco2,aerofixh2o                                        &
82     &   , hydrology, sourceevol                                        &
83     &   , CLFvarying                                                   &
84     &   , strictboundcorrk                                             &                                       
85     &   , ok_slab_ocean                                                &
86     &   , ok_slab_sic                                                  &
87     &   , ok_slab_heat_transp                                         
88
89
90      COMMON/callkeys_i/iaervar,iddist,iradia,startype
91     
92      COMMON/callkeys_r/topdustref,Nmix_co2,dusttau,Fat1AU,stelTbb,     &
93     &                  Tstrat,tplanet,obs_tau_col_tropo,               &
94     &                  obs_tau_col_strato,pres_bottom_tropo,           &
95     &                  pres_top_tropo,pres_bottom_strato,              &
96     &                  pres_top_strato,size_tropo,size_strato,satval,  &
97     &                  CLFfixval,n2mixratio,co2supsat,pceil,albedosnow,&
98     &                  maxicethick,Tsaldiff,tau_relax,cloudlvl,        &
99     &                  icetstep,intheat,flatten,Rmean,J2,MassPlanet
100     
101      logical callrad,corrk,calldifv,UseTurbDiff                        &
102     &   , calladj,co2cond,callsoil                                     &
103     &   , season,diurnal,tlocked,rings_shadow,lwrite                   &
104     &   , callstats,calleofdump                                        &
105     &   , callgasvis,continuum,H2Ocont_simple,graybody                 &
106     &   , strictboundcorrk                                             
107
108      logical enertest
109      logical nonideal
110      logical meanOLR
111      logical specOLR
112      logical kastprof
113      logical newtonian
114      logical check_cpp_match
115      logical force_cpp
116      logical testradtimes
117      logical rayleigh
118      logical stelbbody
119      logical ozone
120      logical nearco2cond
121      logical tracer
122      logical mass_redistrib
123      logical varactive
124      logical varfixed
125      logical radfixed
126      logical sedimentation
127      logical water,watercond,waterrain
128      logical aeroco2,aeroh2o,aeroh2so4,aeroback2lay
129      logical aerofixco2,aerofixh2o
130      logical hydrology
131      logical sourceevol
132      logical CLFvarying
133      logical nosurf
134      logical oblate
135      logical ok_slab_ocean
136      logical ok_slab_sic
137      logical ok_slab_heat_transp
138
139      integer iddist
140      integer iaervar
141      integer iradia
142      integer startype
143
144      real topdustref
145      real Nmix_co2
146      real dusttau
147      real Fat1AU
148      real stelTbb
149      real Tstrat
150      real tplanet
151      real obs_tau_col_tropo
152      real obs_tau_col_strato
153      real pres_bottom_tropo
154      real pres_top_tropo
155      real pres_bottom_strato
156      real pres_top_strato
157      real size_tropo
158      real size_strato
159      real satval
160      real CLFfixval
161      real n2mixratio
162      real co2supsat
163      real pceil
164      real albedosnow
165      real maxicethick
166      real Tsaldiff
167      real tau_relax
168      real cloudlvl
169      real icetstep
170      real intheat
171      real flatten
172      real Rmean
173      real J2
174      real MassPlanet
175!-----------------------------------------------------------------------
176! INCLUDE "comcstfi.h"
177
178      common/comcstfi/pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
179      common/comcstfi/avocado!,molrad,visc
180     
181      real pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
182      real avocado!,molrad,visc
183
184!-----------------------------------------------------------------------
185! INCLUDE planet.h
186
187      COMMON/planet/apoastr,periastr,year_day,peri_day,                 &
188     &       obliquit,nres,                                             &
189     &       z0,lmixmin,emin_turb,coefvis,coefir,                       &
190     &       timeperi,e_elips,p_elips
191
192      real apoastr,periastr,year_day,peri_day,                          &
193     &     obliquit,nres,                                               &
194     &     z0,lmixmin,emin_turb,coefvis,coefir,                         &
195     &       timeperi,e_elips,p_elips
196     
197
198!-----------------------------------------------------------------------
199!#include "control.h"
200
201  ! --------------------------------------------------------------
202  !  Declarations
203  ! --------------------------------------------------------------
204
205  integer nlayer,nlevel,nq
206  integer ilay,ilev,iq,iw,iter
207  real play(nlayermx)     ! Pressure at the middle of the layers [Pa]
208  real zlay(nlayermx)     ! Altitude at middle of the layers [km]
209  real plev(nlayermx+1)   ! Intermediate pressure levels [Pa]
210  real temp(nlayermx)     ! temperature at the middle of the layers [K]
211  real,allocatable :: q(:,:)   ! tracer mixing ratio [kg/kg]
212  real,allocatable :: vmr(:,:) ! tracer mixing ratio [mol/mol]
213  real,allocatable :: qsurf(:)        ! tracer surface budget [kg/kg] ????
214  real psurf,psurf_n,tsurf     
215
216  real emis, albedo
217
218  real muvar(nlayermx+1)
219
220  real dtsw(nlayermx) ! heating rate (K/s) due to SW
221  real dtlw(nlayermx) ! heating rate (K/s) due to LW
222  real fluxsurf_lw   ! incident LW flux to surf (W/m2)
223  real fluxtop_lw    ! outgoing LW flux to space (W/m2)
224  real fluxsurf_sw   ! incident SW flux to surf (W/m2)
225  real fluxabs_sw    ! SW flux absorbed by planet (W/m2)
226  real fluxtop_dn    ! incident top of atmosphere SW flux (W/m2)
227
228  ! not used
229  real reffrad(nlayermx,naerkind)
230  real nueffrad(nlayermx,naerkind)
231  real cloudfrac(nlayermx)
232  real totcloudfrac
233  real tau_col
234
235  real dTstrat
236  real aerosol(nlayermx,naerkind) ! aerosol tau (kg/kg)
237  real OLR_nu(1,L_NSPECTI)
238  real OSR_nu(1,L_NSPECTV)
239  real Eatmtot
240
241  integer ierr
242  logical firstcall,lastcall,global1d
243
244  character*20,allocatable :: nametrac(:)   ! name of the tracer (no need for adv trac common)
245
246  ! --------------
247  ! Initialisation
248  ! --------------
249
250  pi=2.E+0*asin(1.E+0)
251
252  reffrad(:,:)  = 0.0
253  nueffrad(:,:) = 0.0
254  cloudfrac(:)  = 0.0
255  totcloudfrac  = 0.0
256
257
258  nlayer=nlayermx
259  nlevel=nlayer+1
260
261  !! this is defined in comsaison_h
262  ALLOCATE(mu0(1))
263  ALLOCATE(fract(1))
264
265
266
267  !  Load parameters from "run.def"
268  ! -------------------------------
269
270  ! check if 'kcm1d.def' file is around (otherwise reading parameters
271  ! from callphys.def via getin() routine won't work.)
272  open(90,file='kcm1d.def',status='old',form='formatted',&
273       iostat=ierr)
274  if (ierr.ne.0) then
275     write(*,*) 'Cannot find required file "kcm1d.def"'
276     write(*,*) '  (which should contain some input parameters'
277     write(*,*) '   along with the following line:'
278     write(*,*) '   INCLUDEDEF=callphys.def'
279     write(*,*) '   )'
280     write(*,*) ' ... might as well stop here ...'
281     stop
282  else
283     close(90)
284  endif
285
286! now, run.def is needed anyway. so we create a dummy temporary one
287! for ioipsl to work. if a run.def is already here, stop the
288! program and ask the user to do a bit of cleaning
289  open(90,file='run.def',status='old',form='formatted',& 
290       iostat=ierr)
291  if (ierr.eq.0) then
292     close(90)
293     write(*,*) 'There is already a run.def file.'
294     write(*,*) 'This is not compatible with 1D runs.'
295     write(*,*) 'Please remove the file and restart the run.'
296     write(*,*) 'Runtime parameters are supposed to be in kcm1d.def'
297     stop
298  else
299     call system('touch run.def')
300     call system("echo 'INCLUDEDEF=callphys.def' >> run.def")
301     call system("echo 'INCLUDEDEF=kcm1d.def' >> run.def")
302  endif
303
304  global1d = .false. ! default value
305  call getin("global1d",global1d)
306  if(.not.global1d)then
307     print*,'Error, kcm1d must have global1d=.true. in kcm1d.def!'
308     stop
309  end if
310
311  psurf_n=100000. ! default value for psurf
312  print*,'Dry surface pressure (Pa)?'
313  call getin("psurf",psurf_n)
314  write(*,*) " psurf = ",psurf_n
315
316! OK. now that run.def has been read once -- any variable is in memory.
317! so we can dump the dummy run.def
318  call system("rm -rf run.def")
319
320  tsurf=300.0 ! default value for tsurf
321  print*,'Surface temperature (K)?'
322  call getin("tref",tsurf)
323  write(*,*) " tsurf = ",tsurf
324
325  g=10.0 ! default value for g
326  print*,'Gravity ?'
327  call getin("g",g)
328  write(*,*) " g = ",g
329
330  periastr = 1.0
331  apoastr  = 1.0
332  print*,'Periastron (AU)?'
333  call getin("periastr",periastr)
334  write(*,*) "periastron = ",periastr
335  dist_star = periastr 
336  fract     = 0.5
337  print*,'Apoastron (AU)?'
338  call getin("apoastr",apoastr)
339  write(*,*) "apoastron = ",apoastr
340
341  albedo=0.2 ! default value for albedo
342  print*,'Albedo of bare ground?'
343  call getin("albedo",albedo)
344  write(*,*) " albedo = ",albedo
345
346  emis=1.0 ! default value for emissivity
347  PRINT *,'Emissivity of bare ground ?'
348  call getin("emis",emis)
349  write(*,*) " emis = ",emis
350
351  pceil=100.0 ! Pascals
352  PRINT *,'Ceiling pressure (Pa) ?'
353  call getin("pceil",pceil)
354  write(*,*) " pceil = ", pceil
355
356  mugaz=0.
357  cpp=0.
358
359  check_cpp_match = .false.
360  call getin("check_cpp_match",check_cpp_match)
361  if (check_cpp_match) then
362     print*,"In 1D modeling, check_cpp_match is supposed to be F"
363     print*,"Please correct callphys.def"
364     stop
365  endif
366
367  call su_gases
368  call calc_cpp_mugaz
369
370  call inifis(1,llm,0,86400.0,1.0,0.0,0.0,1.0,rad,g,r,cpp)
371
372  ! Tracer initialisation
373  ! ---------------------
374  if (tracer) then
375     ! load tracer names from file 'traceur.def'
376     open(90,file='traceur.def',status='old',form='formatted',&
377          iostat=ierr)
378     if (ierr.eq.0) then
379        write(*,*) "kcm1d: Reading file traceur.def"
380        ! read number of tracers:
381        read(90,*,iostat=ierr) nq
382        if (ierr.ne.0) then
383           write(*,*) "kcm1d: error reading number of tracers"
384           write(*,*) "   (first line of traceur.def) "
385           stop
386        endif
387        nqtot=nq
388        ! allocate arrays which depend on number of tracers
389        allocate(nametrac(nq))
390        allocate(q(nlayermx,nq))
391        allocate(vmr(nlayermx,nq))
392        allocate(qsurf(nq))
393
394        do iq=1,nq
395           ! minimal version, just read in the tracer names, 1 per line
396           read(90,*,iostat=ierr) nametrac(iq)
397           if (ierr.ne.0) then
398              write(*,*) 'kcm1d: error reading tracer names...'
399              stop
400           endif
401        enddo !of do iq=1,nq
402     endif
403
404     call initracer(1,nq,nametrac)
405
406  endif
407
408
409  do iq=1,nq
410     do ilay=1,nlayer
411        q(ilay,iq) = 0.
412     enddo
413  enddo
414
415  do iq=1,nq
416     qsurf(iq) = 0.
417  enddo
418
419  firstcall = .true.
420  lastcall  = .false.
421
422  iter    = 1
423  Tstrat  = 150.0
424  dTstrat = 1000.0
425
426  ! ---------
427  ! Run model
428  ! ---------
429  !do
430     psurf = psurf_n
431
432     !    Create vertical profiles
433     call kcmprof_fn(psurf,qsurf(1),tsurf,    &
434          tstrat,play,plev,zlay,temp,q(:,1),muvar(1))
435
436     !    Run radiative transfer
437     call callcorrk(1,nlayer,q,nq,qsurf,      &
438          albedo,emis,mu0,plev,play,temp,                    &
439          tsurf,fract,dist_star,aerosol,muvar,         &
440          dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
441          fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,reffrad,nueffrad,tau_col,  &
442          cloudfrac,totcloudfrac,.false.,firstcall,lastcall)
443
444     !    Iterate stratospheric temperature
445     print*,'Tstrat = ',Tstrat
446     dTstrat = Tstrat
447     !Tstrat  = Tsurf*(0.2786*(psurf/100000.)**(-1.123))**0.25
448     ! skin temperature (gray approx.) using analytic pure H2 expression
449     !Tstrat  = (fluxabs_sw/(2*sigma))**0.25 ! skin temperature (gray approx.)
450     Tstrat  = (fluxtop_lw/(2*sigma))**0.25 ! skin temperature (gray approx.)
451     dTstrat = dTstrat-Tstrat
452
453     !if(abs(dTstrat).lt.1.0)then
454     !   print*,'dTstrat = ',dTstrat
455     !   exit
456     !endif
457
458     !iter=iter+1
459     !if(iter.eq.100)then
460     !   print*,'Stratosphere failed to converge after'
461     !   print*,'100 iteration, aborting run.'
462     !   call abort
463     !endif
464
465  !end do
466
467  ! Run radiative transfer one last time to get OLR,OSR
468  firstcall=.false.
469  lastcall=.true.
470  call callcorrk(1,nlayer,q,nq,qsurf,      &
471       albedo,emis,mu0,plev,play,temp,                    &
472       tsurf,fract,dist_star,aerosol,muvar,         &
473       dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
474       fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,             &
475       reffrad,nueffrad,tau_col,  &
476       cloudfrac,totcloudfrac,.false.,firstcall,lastcall)
477
478
479  ! Calculate total atmospheric energy
480  Eatmtot=0.0
481  !  call calcenergy_kcm(tsurf,temp,play,plev,qsurf,&
482  !     q(:,1),muvar,Eatmtot)
483
484  ! ------------------------
485  ! Save data to ascii files
486  ! ------------------------
487
488  print*,'Saving profiles...'
489  open(115,file='profpres.out',form='formatted')
490  open(116,file='proftemp.out',form='formatted')
491  open(117,file='profztab.out',form='formatted')
492  open(118,file='profqvar.out',form='formatted')
493  open(119,file='profvvar.out',form='formatted')
494
495  write(115,*) psurf
496  write(116,*) tsurf
497  write(117,*) 0.0
498  write(118,*) qsurf(1)
499  write(119,*) qsurf(1)*(muvar(1)/mH2O)
500  do ilay=1,nlayer
501     vmr(ilay,1) = q(ilay,1)*(muvar(ilay+1)/mH2O)
502     write(115,*) play(ilay)
503     write(116,*) temp(ilay)
504     write(117,*) zlay(ilay)
505     write(118,*) q(ilay,1)
506     write(119,*) vmr(ilay,1)
507  enddo
508  close(115)
509  close(116)
510  close(117)
511  close(118)
512  close(119)
513
514  print*, tsurf,psurf,fluxtop_dn,fluxabs_sw,fluxtop_lw 
515
516  print*,'Saving scalars...'
517  open(116,file='surf_vals.out')
518  write(116,*) tsurf,psurf,fluxtop_dn,         &
519       fluxabs_sw,fluxtop_lw 
520  close(116)
521  open(111,file='ene_vals.out')
522  write(111,*) tsurf,psurf,Eatmtot,Tstrat
523  close(111)
524
525  print*,'Saving spectra...'
526  open(117,file='OLRnu.out')
527  do iw=1,L_NSPECTI
528     write(117,*) OLR_nu(1,iw)
529  enddo
530  close(117)
531
532  open(127,file='OSRnu.out')
533  do iw=1,L_NSPECTV
534     write(127,*) OSR_nu(1,iw)
535  enddo
536  close(127) 
537
538end program kcm1d
Note: See TracBrowser for help on using the repository browser.