source: trunk/LMDZ.GENERIC/libf/phystd/dyn1d/kcm1d.F90 @ 3574

program kcm1d

  use infotrac, only: nqtot
  use radinc_h,      only: NAERKIND
  use radcommon_h,   only: L_NSPECTI, L_NSPECTV, sigma, glat
  use watercommon_h, only: mH2O
  use ioipsl_getincom, only: getin
  use comsaison_h, only: mu0, fract, dist_star
  use planete_mod
  use callkeys_mod, only: pceil, tstrat, tracer, global1d, &
                          varspec, varspec_data, nvarlayer
  use inifis_mod, only: inifis
  use aerosol_mod, only: iniaerosol
  use callcorrk_mod, only: callcorrk
  use comcstfi_mod
  use mod_grid_phy_lmdz, only : regular_lonlat
  use regular_lonlat_mod, only: init_regular_lonlat
  use physics_distribution_mod, only: init_physics_distribution
  use regular_lonlat_mod, only: init_regular_lonlat
  use mod_interface_dyn_phys, only: init_interface_dyn_phys
  use geometry_mod, only: init_geometry
  use dimphy, only : init_dimphy
  use gases_h, only: ngasmx
  use version_info_mod, only: print_version_info

  implicit none

  !==================================================================
  !     
  !     Purpose
  !     -------
  !     Run the universal model radiative transfer once in a 1D column.
  !     Useful for climate evolution studies etc.
  !     
  !     It can be compiled with a command like (e.g. for 25 layers):
  !                                  "makegcm -p std -d 25 kcm1d"
  !     It requires the files "callphys.def", "gases.def"
  !     "traceur.def", and "run.def" with a line "INCLUDEDEF=callphys.def"
  !
  !     Authors
  !     -------
  !     - R. Wordsworth
  !     - updated by M. Turbet (June 2017)
  !
  !==================================================================

#include "dimensions.h"
!#include "dimphys.h"

  ! --------------------------------------------------------------
  !  Declarations
  ! --------------------------------------------------------------

  integer nlayer,nlevel,nq
  integer ilay,ilev,iq,iw,iter
  real play(llm)     ! Pressure at the middle of the layers [Pa]
  real zlay(llm)     ! Altitude at middle of the layers [km]
  real plev(llm+1)   ! Intermediate pressure levels [Pa]
  real temp(llm)     ! temperature at the middle of the layers [K]
  real,allocatable :: q(:,:)   ! tracer mixing ratio [kg/kg]
  real,allocatable :: vmr(:,:) ! tracer mixing ratio [mol/mol]
  real,allocatable :: qsurf(:)        ! tracer surface budget [kg/kg] ????
  real psurf,psurf_n,tsurf(1)      

  real emis(1), albedo(1), albedo_equivalent(1)
  real albedo_wv(1,L_NSPECTV)

  real muvar(llm+1)

  real dtsw(llm) ! heating rate (K/s) due to SW
  real dtlw(llm) ! heating rate (K/s) due to LW
  real fluxsurf_lw(1)   ! incident LW flux to surf (W/m2)
  real fluxtop_lw(1)    ! outgoing LW flux to space (W/m2)
  real fluxsurf_sw(1)   ! incident SW flux to surf (W/m2)
  real fluxsurfabs_sw(1) ! absorbed SW flux by the surf (W/m2)
  real fluxabs_sw(1)    ! SW flux absorbed by planet (W/m2)
  real fluxtop_dn(1)    ! incident top of atmosphere SW flux (W/m2)

  ! not used
  real cloudfrac(1,llm)
  real totcloudfrac(1)
  real tau_col(1)

  real dTstrat
  real,allocatable :: aerosol(:,:) ! aerosol tau (kg/kg)
  real OLR_nu(1,L_NSPECTI)
  real OSR_nu(1,L_NSPECTV)
  real GSR_nu(1,L_NSPECTV)
  real int_dtaui(1,llm,L_NSPECTI)
  real int_dtauv(1,llm,L_NSPECTV)
  real Eatmtot
  
  ! For fixed variable molar mass
  real, dimension(:),allocatable,save :: p_var
  real, dimension(:),allocatable,save :: mu_var
  real, dimension(:,:),allocatable,save :: frac_var

  integer ierr
  logical firstcall,lastcall

  character*20,allocatable :: nametrac(:)   ! name of the tracer (no need for adv trac common)
  
  real :: latitude(1), longitude(1), cell_area(1), phisfi(1)
  
  !     added by JVO and YJ to read modern traceur.def
  character(len=500) :: line ! to store a line of text
  LOGICAL :: moderntracdef=.false. ! JVO, YJ : modern traceur.def

  character(len=100) :: dt_file
  integer :: ios
  integer :: k

  character(100) :: arg ! To read command-line arguments
  
  ! --------------
  ! Initialisation
  ! --------------
  if (command_argument_count() > 0) then ! Get the number of command-line arguments
      call get_command_argument(1,arg) ! Read the argument given to the program
      select case (trim(adjustl(arg)))
          case('version')
              call print_version_info()
              stop
          case default
              error stop 'The argument given to the program is unknown!'
      end select
  endif

  pi=2.E+0*asin(1.E+0)

  cloudfrac(1,:)  = 0.0
  totcloudfrac(1)  = 0.0


  nlayer=llm
  nlevel=nlayer+1

  !! this is defined in comsaison_h
  ALLOCATE(mu0(1))
  ALLOCATE(fract(1))
  ALLOCATE(glat(1))

! Initialize fixed variable mu
!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~   

  if(varspec) then
    IF (.NOT.ALLOCATED(p_var))  ALLOCATE(p_var(nvarlayer))
    IF (.NOT.ALLOCATED(mu_var))  ALLOCATE(mu_var(nvarlayer))
    IF (.NOT.ALLOCATED(frac_var))  ALLOCATE(frac_var(nvarlayer,ngasmx))
    p_var = 0.0
    mu_var = 0.0
    frac_var = 0.0
    dt_file=TRIM(varspec_data)
    open(34,file=dt_file,form='formatted',status='old',iostat=ios)
    if (ios.ne.0) then        ! file not found
      write(*,*) 'Error from varspec'
      write(*,*) 'Data file ',trim(varspec_data),' not found.'
      write(*,*) 'Check that the data is in your run repository.'
      call abort_physic !a verifier
    else
      do k=1,nvarlayer
        read(34,*) p_var(k), mu_var(k),frac_var(k,1:ngasmx)
        !The order of columns in frac_var must correspond to order of molecules gases.def 
        !The format of your file must be:
        ! pressure(k) molar_mass(k), molar_fraction_of_gas_1(k), molar_fraction_of_gas_2(k),..., molar_fraction_of_gas_ngasmx(k)
      enddo
    endif
    close(34)
  endif

  !  Load parameters from "run.def" 
  ! -------------------------------

  ! check if 'kcm1d.def' file is around (otherwise reading parameters
  ! from callphys.def via getin() routine won't work.)
  open(90,file='kcm1d.def',status='old',form='formatted',&
       iostat=ierr)
  if (ierr.ne.0) then
     write(*,*) 'Cannot find required file "kcm1d.def"'
     write(*,*) '  (which should contain some input parameters'
     write(*,*) '   along with the following line:'
     write(*,*) '   INCLUDEDEF=callphys.def'
     write(*,*) '   )'
     write(*,*) ' ... might as well stop here ...'
     stop
  else
     close(90)
  endif
  
! now, run.def is needed anyway. so we create a dummy temporary one
! for ioipsl to work. if a run.def is already here, stop the
! program and ask the user to do a bit of cleaning
  open(90,file='run.def',status='old',form='formatted',& 
       iostat=ierr)
  if (ierr.eq.0) then
     close(90)
     write(*,*) 'There is already a run.def file.'
     write(*,*) 'This is not compatible with 1D runs.'
     write(*,*) 'Please remove the file and restart the run.'
     write(*,*) 'Runtime parameters are supposed to be in kcm1d.def'
     stop
  else
     call system('touch run.def')
     call system("echo 'INCLUDEDEF=callphys.def' >> run.def")
     call system("echo 'INCLUDEDEF=kcm1d.def' >> run.def")
  endif
    
! check if we are going to run with or without tracers
  write(*,*) "Run with or without tracer transport ?"
  tracer=.false. ! default value
  call getin("tracer",tracer)
  write(*,*) " tracer = ",tracer



  ! Tracer initialisation
  ! ---------------------
  if (tracer) then
     ! load tracer names from file 'traceur.def'
     open(90,file='traceur.def',status='old',form='formatted',&
          iostat=ierr)
     if (ierr.eq.0) then
        write(*,*) "kcm1d: Reading file traceur.def"
        ! read number of tracers:
        !! - Modif. by JVO and YJ for modern planetary traceur.def ---------------
        READ(90,'(A)') line
        IF (trim(line).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
           READ(line,*,iostat=ierr) nq ! Try standard traceur.def
        ELSE
           moderntracdef = .true.
           DO
             READ(90,'(A)',iostat=ierr) line
             IF (ierr.eq.0) THEN
               IF (index(line,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
                 READ(line,*,iostat=ierr) nq
                 EXIT
               ENDIF
             ELSE ! If pb, or if reached EOF without having found nbtr
                write(*,*) "rcm1d: error reading number of tracers"
                write(*,*) "   (first line of traceur.def) "
                stop
             ENDIF
           ENDDO
        ENDIF ! if modern or standard traceur.def
        if (ierr.ne.0) then
           write(*,*) "kcm1d: error reading number of tracers"
           write(*,*) "   (first line of traceur.def) "
           stop
        endif
        nqtot=nq
        ! allocate arrays which depend on number of tracers
        allocate(nametrac(nq))
        allocate(q(nlayer,nq))
        allocate(vmr(nlayer,nq))
        allocate(qsurf(nq))

        do iq=1,nq
           ! minimal version, just read in the tracer names, 1 per line
           read(90,*,iostat=ierr) nametrac(iq)
           write(*,*) 'tracer here is', nametrac(iq)
           if (ierr.ne.0) then
              write(*,*) 'kcm1d: error reading tracer names...'
              stop
           endif
        enddo !of do iq=1,nq
        close(90)
     endif

     call initracer(1,nq)

  endif
  




  psurf_n=100000. ! default value for psurf
  print*,'Dry surface pressure (Pa)?'
  call getin("psurf",psurf_n)
  write(*,*) " psurf = ",psurf_n

! OK. now that run.def has been read once -- any variable is in memory.
! so we can dump the dummy run.def
  call system("rm -rf run.def")

  tsurf=300.0 ! default value for tsurf
  print*,'Surface temperature (K)?'
  call getin("tref",tsurf)
  write(*,*) " tsurf = ",tsurf

  g=10.0 ! default value for g
  print*,'Gravity ?'
  call getin("g",g)
  write(*,*) " g = ",g
  glat(1)=g

  periastr = 1.0
  apoastr  = 1.0
  print*,'Periastron (AU)?'
  call getin("periastr",periastr)
  write(*,*) "periastron = ",periastr
  dist_star = periastr 
  fract     = 0.5
  print*,'Apoastron (AU)?'
  call getin("apoastr",apoastr)
  write(*,*) "apoastron = ",apoastr

  albedo=0.2 ! default value for albedo
  print*,'Albedo of bare ground?'
  call getin("albedo",albedo)
  write(*,*) " albedo = ",albedo
  do iw=1,L_NSPECTV
     albedo_wv(1,iw)=albedo(1)
  enddo

  emis=1.0 ! default value for emissivity
  PRINT *,'Emissivity of bare ground ?'
  call getin("emis",emis)
  write(*,*) " emis = ",emis

  pceil=100.0 ! Pascals
  PRINT *,'Ceiling pressure (Pa) ?'
  call getin("pceil",pceil)
  write(*,*) " pceil = ", pceil

!!! GEOGRAPHICAL INITIALIZATIONS
     !!! AREA. useless in 1D
  cell_area(1)=1.E+0 !JL+EM to have access to the area in the diagfi.nc files.
     !!! GEOPOTENTIAL. useless in 1D because control by surface pressure
  phisfi(1)=0.E+0
     !!! LATITUDE. can be set. 
  latitude=0 ! default value for latitude
  PRINT *,'latitude (in degrees) ?'
  call getin("latitude",latitude)
  write(*,*) " latitude = ",latitude
  latitude=latitude*pi/180.E+0
     !!! LONGITUDE. useless in 1D.
  longitude=0.E+0
  longitude=longitude*pi/180.E+0

  rad=6400000.
  !rad = -99999.
  !PRINT *,'PLANETARY RADIUS in m ?'
  !call getin("rad",rad)
  ! Planetary  radius is needed to compute shadow of the rings
  !IF (rad.eq.-99999. .and. rings_shadow .eqv. .true.) THEN
  !   PRINT *,"STOP. I NEED rad IN RCM1D.DEF."
  !   STOP
  !ELSE
  !   PRINT *,"--> rad = ",rad
  !ENDIF



  call init_physics_distribution(regular_lonlat,4,1,1,1,nlayer,1)
  call init_interface_dyn_phys
  CALL init_regular_lonlat(1,1,longitude,latitude,(/0.,0./),(/0.,0./))
  call init_geometry(1,longitude,latitude,(/0.,0.,0.,0./),(/0.,0.,0.,0./),cell_area)
  call init_dimphy(1,nlayer) ! Initialize dimphy module

!!! CALL INIFIS
!!! - read callphys.def
!!! - calculate sine and cosine of longitude and latitude
!!! - calculate mugaz and cp 
!!! NB: some operations are being done dummily in inifis in 1D
!!! - physical constants: nevermind, things are done allright below
!!! - physical frequency: nevermind, in inifis this is a simple print

  cpp=1.d-7 !JL because we divide by cpp in inifis, there may be a more elegant solution
  CALL inifis(1,llm,nq,0,86400.0,1,1.0,latitude,longitude,cell_area,rad,g,r,cpp)

  if(.not.global1d)then
     print*,'Error, kcm1d must have global1d=.true. in kcm1d.def!'
     stop
  end if

  !write(*,*) 'BASE 1'
  !write(*,*) 1,llm,nq,0,86400.0,1,1.0,latitude,longitude,cell_area,rad,g,r,cpp

  ! initialise naerkind (from callphys.def) and allocate aerosol(:,:)
  naerkind=0 !default
  call getin("naerkind",naerkind)
  allocate(aerosol(llm,naerkind))
  aerosol(:,:)=0

  do iq=1,nq
     do ilay=1,nlayer
        q(ilay,iq) = 0.
     enddo
  enddo
  
  do iq=1,nq
     qsurf(iq) = 0.
  enddo

  firstcall = .true.
  lastcall  = .false.

  iter    = 1
  Tstrat  = 150.0
  dTstrat = 1000.0

  ! ---------
  ! Run model
  ! ---------
  !do
  psurf = psurf_n

     !    Create vertical profiles
  call kcmprof_fn(nlayer,psurf,qsurf(1),tsurf,    &
       tstrat,play,plev,zlay,temp,q(:,1),muvar(1))
         

  !write(*,*) 'BASE 2'
  !write(*,*) nlayer,psurf,qsurf(1),tsurf,    &
  !           tstrat
  !write(*,*) play
  !write(*,*) plev
  !write(*,*) zlay
  !write(*,*) temp
  !write(*,*) q(:,1),muvar(1)

         
               
     call iniaerosol
     

     

     !    Run radiative transfer
     call callcorrk(1,nlayer,q,nq,qsurf,                  &
          albedo_wv,albedo_equivalent,                    &
          emis,mu0,plev,play,temp,                        &
          tsurf,fract,dist_star,aerosol,muvar,            &
          dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,              & 
          fluxsurfabs_sw,fluxtop_lw,                      &
          fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,GSR_nu,     &
          int_dtaui,int_dtauv,                            &
          tau_col,cloudfrac,totcloudfrac,                 &
          .false.,p_var,frac_var,firstcall,lastcall)

     !write(*,*) 'BASE 3'
     !write(*,*) 1,nlayer,'A',q,'B',nq,'C',qsurf,                      &
     !     albedo_wv,'D',albedo_equivalent,'E',                    &
     !     emis,'F',mu0,'G',plev,'H',play,'I',temp,'J',                        &
     !     tsurf,'K',fract,'L',dist_star,'M',aerosol,'N',muvar,'O',            &
     !     dtlw,'P',dtsw,'Q',fluxsurf_lw,'R',fluxsurf_sw,'S',              & 
     !     fluxsurfabs_sw,'T',fluxtop_lw,'U',                      &
     !     fluxabs_sw,'V',fluxtop_dn,'W',OLR_nu,'X',OSR_nu,'Y',tau_col,'Z',    &
     !     cloudfrac,'A1',totcloudfrac,'A2',.false.,'A3',firstcall,'A4',lastcall







     !    Iterate stratospheric temperature
     print*,'Tstrat = ',Tstrat
     dTstrat = Tstrat
     !Tstrat  = Tsurf*(0.2786*(psurf/100000.)**(-1.123))**0.25 
     ! skin temperature (gray approx.) using analytic pure H2 expression
     !Tstrat  = (fluxabs_sw/(2*sigma))**0.25 ! skin temperature (gray approx.)
     Tstrat  = (fluxtop_lw(1)/(2*sigma))**0.25 ! skin temperature (gray approx.)
     dTstrat = dTstrat-Tstrat

     !if(abs(dTstrat).lt.1.0)then
     !   print*,'dTstrat = ',dTstrat
     !   exit
     !endif

     !iter=iter+1
     !if(iter.eq.100)then
     !   print*,'Stratosphere failed to converge after'
     !   print*,'100 iteration, aborting run.'
     !   call abort
     !endif

  !end do


  ! Run radiative transfer one last time to get OLR,OSR
  firstcall=.false.
  lastcall=.true.
  call callcorrk(1,nlayer,q,nq,qsurf,                          &
       albedo_wv,albedo_equivalent,emis,mu0,plev,play,temp,    &
       tsurf,fract,dist_star,aerosol,muvar,                    &
       dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxsurfabs_sw,       &
       fluxtop_lw, fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,GSR_nu, &
       int_dtaui,int_dtauv,                                    &
       tau_col,cloudfrac,totcloudfrac,                         &
       .false.,p_var,frac_var,firstcall,lastcall)


     !write(*,*) 'BASE 4'
     !write(*,*) 1,nlayer,'A',q,'B',nq,'C',qsurf,                      &
     !     albedo_wv,'D',albedo_equivalent,'E',                    &
     !     emis,'F',mu0,'G',plev,'H',play,'I',temp,'J',                        &
     !     tsurf,'K',fract,'L',dist_star,'M',aerosol,'N',muvar,'O',            &
     !     dtlw,'P',dtsw,'Q',fluxsurf_lw,'R',fluxsurf_sw,'S',              & 
     !     fluxsurfabs_sw,'T',fluxtop_lw,'U',                      &
     !     fluxabs_sw,'V',fluxtop_dn,'W',OLR_nu,'X',OSR_nu,'Y',tau_col,'Z',    &
     !     cloudfrac,'A1',totcloudfrac,'A2',.false.,'A3',firstcall,'A4',lastcall







  ! Calculate total atmospheric energy
  Eatmtot=0.0
  !  call calcenergy_kcm(nlayer,tsurf,temp,play,plev,qsurf,&
  !     q(:,1),muvar,Eatmtot)

  ! ------------------------
  ! Save data to ascii files
  ! ------------------------

  print*,'Saving profiles...'
  open(115,file='profpres.out',form='formatted')
  open(116,file='proftemp.out',form='formatted')
  open(117,file='profztab.out',form='formatted')
  open(118,file='profqvar.out',form='formatted')
  open(119,file='profvvar.out',form='formatted')

  write(115,*) psurf
  write(116,*) tsurf
  write(117,*) 0.0
  write(118,*) qsurf(1)
  write(119,*) qsurf(1)*(muvar(1)/mH2O)
  do ilay=1,nlayer
     vmr(ilay,1) = q(ilay,1)*(muvar(ilay+1)/mH2O)
     write(115,*) play(ilay)
     write(116,*) temp(ilay)
     write(117,*) zlay(ilay)
     write(118,*) q(ilay,1)
     write(119,*) vmr(ilay,1)
  enddo
  close(115)
  close(116)
  close(117)
  close(118)
  close(119)

  print*, tsurf,psurf,fluxtop_dn,fluxabs_sw,fluxtop_lw 

  print*,'Saving scalars...'
  open(116,file='surf_vals.out')
  write(116,*) tsurf,psurf,fluxtop_dn,         &
       fluxabs_sw,fluxtop_lw 
  close(116)
  open(111,file='ene_vals.out')
  write(111,*) tsurf,psurf,Eatmtot,Tstrat
  close(111)

  print*,'Saving spectra...'
  open(117,file='OLRnu.out')
  do iw=1,L_NSPECTI
     write(117,*) OLR_nu(1,iw)
  enddo
  close(117)

  open(127,file='OSRnu.out')
  do iw=1,L_NSPECTV
     write(127,*) OSR_nu(1,iw)
  enddo
  close(127)  

end program kcm1d