Changeset 1714

Timestamp:

Jun 7, 2017, 4:05:31 PM (8 years ago)

Author:

mturbet

Message:

kcm1d like a phoenix rises from the ashes, part I

Location:

trunk/LMDZ.GENERIC

Files:

• README (modified) (1 diff)
• libf/phystd/callcorrk.F90 (modified) (2 diffs)
• libf/phystd/dyn1d/gradients_kcm.F90 (modified) (1 diff)
• libf/phystd/dyn1d/kcm1d.F90 (modified) (15 diffs)
• libf/phystd/dyn1d/kcmprof_fn.F90 (modified) (1 diff)
• libf/phystd/suaer_corrk.F90 (modified) (2 diffs)

trunk/LMDZ.GENERIC/README

@@ -1345 +1345 @@
 Added surfalbedo and surfemis keywords to be used with startphy_file = .false.
 Also made default values in phyetat0_mod unambiguously float
+
+== 07/06/2017 == MT
+Resurrection of kcm1d, part I

trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90

@@ -323 +323 @@
      end do !iaer=1,naerkind.
 
-
       ! How much light do we get ?
       do nw=1,L_NSPECTV
@@ -575 +574 @@
 
       if(kastprof)then
-
+      
+         if(.not.global1d)then ! garde-fou/safeguard added by MT (to be removed in the future)
+            write(*,*) 'You have to fix mu0, '
+            write(*,*) 'the cosinus of the solar angle'
+            stop
+         endif
+         
          ! Initial values equivalent to mugaz.
          DO l=1,nlayer

trunk/LMDZ.GENERIC/libf/phystd/dyn1d/gradients_kcm.F90

@@ -116 +116 @@
         cp_v = 2.226d3
      endif
-
+     
      dTdp  = 1/(rho_n*cp_n+rho_v*cp_v)
      dPndp = 1/(1d0+m_n/m_v*rho_v/rho_n)

trunk/LMDZ.GENERIC/libf/phystd/dyn1d/kcm1d.F90

@@ -3 +3 @@
   use infotrac, only: nqtot
   use radinc_h,      only: NAERKIND
-  use radcommon_h,   only: L_NSPECTI, L_NSPECTV, sigma
+  use radcommon_h,   only: L_NSPECTI, L_NSPECTV, sigma, glat
   use watercommon_h, only: mH2O
   use ioipsl_getincom, only: getin
@@ -11 +11 @@
   use inifis_mod, only: inifis
   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
+
   implicit none
 
@@ -27 +35 @@
   !     Authors
   !     -------
-  !     R. Wordsworth
+  !     - R. Wordsworth
+  !     - updated by M. Turbet (June 2017)
   !
   !==================================================================
@@ -49 +58 @@
   real psurf,psurf_n,tsurf      
 
-  real emis, albedo
+  real emis, albedo, albedo_equivalent
+  real albedo_wv(1,L_NSPECTV)
 
   real muvar(llm+1)
@@ -58 +68 @@
   real fluxtop_lw    ! outgoing LW flux to space (W/m2)
   real fluxsurf_sw   ! incident SW flux to surf (W/m2)
+  real fluxsurfabs_sw ! absorbed SW flux by the surf (W/m2)
   real fluxabs_sw    ! SW flux absorbed by planet (W/m2)
   real fluxtop_dn    ! incident top of atmosphere SW flux (W/m2)
 
   ! not used
-  real reffrad(llm,naerkind)
-  real nueffrad(llm,naerkind)
   real cloudfrac(llm)
   real totcloudfrac
@@ -78 +87 @@
 
   character*20,allocatable :: nametrac(:)   ! name of the tracer (no need for adv trac common)
-
+  
+  real :: latitude(1), longitude(1), cell_area(1), phisfi(1)
+  
   ! --------------
   ! Initialisation
@@ -85 +96 @@
   pi=2.E+0*asin(1.E+0)
 
-  reffrad(:,:)  = 0.0
-  nueffrad(:,:) = 0.0
   cloudfrac(:)  = 0.0
   totcloudfrac  = 0.0
@@ -97 +106 @@
   ALLOCATE(mu0(1))
   ALLOCATE(fract(1))
-
-
+  ALLOCATE(glat(1))
 
   !  Load parameters from "run.def" 
@@ -118 +126 @@
      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
@@ -136 +144 @@
      call system("echo 'INCLUDEDEF=kcm1d.def' >> run.def")
   endif
-
-  global1d = .false. ! default value
-  call getin("global1d",global1d)
-  if(.not.global1d)then
-     print*,'Error, kcm1d must have global1d=.true. in kcm1d.def!'
-     stop
-  end if
-
-  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
-
-  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
-
-  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
-
-  mugaz=0.
-  cpp=0.
-
-  check_cpp_match = .false.
-  call getin("check_cpp_match",check_cpp_match)
-  if (check_cpp_match) then
-     print*,"In 1D modeling, check_cpp_match is supposed to be F"
-     print*,"Please correct callphys.def"
-     stop
-  endif
-
-  call su_gases
-  call calc_cpp_mugaz
-
-  call inifis(1,llm,0,1,86400.0,1,1.0,&
-              (/0.0/),(/0.0/),(/1.0/),&
-              rad,g,r,cpp)
+    
+! 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
@@ -232 +178 @@
            ! 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...'
@@ -242 +189 @@
 
   endif
-
+  
+
+
+
+  global1d = .false. ! default value
+  call getin("global1d",global1d)
+  write(*,*) " global1d = ",global1d
+  if(.not.global1d)then
+     print*,'Error, kcm1d must have global1d=.true. in kcm1d.def!'
+     stop
+  end if
+
+  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
+  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
+
+  check_cpp_match = .false.
+  call getin("check_cpp_match",check_cpp_match)
+  if (check_cpp_match) then
+     print*,"In 1D modeling, check_cpp_match is supposed to be F"
+     print*,"Please correct callphys.def"
+     stop
+  endif
+
+
+!!! 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)
+
+  !write(*,*) 'BASE 1'
+  !write(*,*) 1,llm,nq,0,86400.0,1,1.0,latitude,longitude,cell_area,rad,g,r,cpp
 
   do iq=1,nq
@@ -249 +313 @@
      enddo
   enddo
-
+  
   do iq=1,nq
      qsurf(iq) = 0.
@@ -265 +329 @@
   ! ---------
   !do
-     psurf = psurf_n
+  psurf = psurf_n
 
      !    Create vertical profiles
-     call kcmprof_fn(nlayer,psurf,qsurf(1),tsurf,    &
-          tstrat,play,plev,zlay,temp,q(:,1),muvar(1))
+  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,emis,mu0,plev,play,temp,                    &
-          tsurf,fract,dist_star,aerosol,muvar,         &
-          dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
-          fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,reffrad,nueffrad,tau_col,  &
+     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,tau_col,    &
           cloudfrac,totcloudfrac,.false.,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
@@ -302 +401 @@
   !end do
 
+
   ! Run radiative transfer one last time to get OLR,OSR
   firstcall=.false.
   lastcall=.true.
-  call callcorrk(1,nlayer,q,nq,qsurf,      &
-       albedo,emis,mu0,plev,play,temp,                    &
-       tsurf,fract,dist_star,aerosol,muvar,         &
-       dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
-       fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,             &
-       reffrad,nueffrad,tau_col,  &
+  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,        &
+       tau_col,                                                &
        cloudfrac,totcloudfrac,.false.,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
+
+
+
+
+
 
 

trunk/LMDZ.GENERIC/libf/phystd/dyn1d/kcmprof_fn.F90

@@ -183 +183 @@
   ! log pressure grid spacing (constant)
   dlogp = -(log(Pn(1)+Pv(1))-log(ptop))/(nlay-1)
-
+  
   call gradients_kcm(profil_flag(1),rho_v(1),rho_n(1),Tsurf,dTdp,dPvdp,dPndp)
   if(verbose)then

trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90

@@ -100 +100 @@
 !$OMP THREADPRIVATE(file_id)
 !---- Please indicate the names of the optical property files below
-!     Please also choose the reference wavelengths of each aerosol
+!     Please also choose the reference wavelengths of each aerosol      
 
       if (noaero) then
@@ -230 +230 @@
 
 !$OMP MASTER
+
             INQUIRE(FILE=TRIM(datadir)//'/'//TRIM(aerdir)//&
                     '/'//TRIM(file_id(iaer,idomain)),&