Ignore:
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/libf/phystd/dyn1d
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r1403 r1714  
    116116        cp_v = 2.226d3
    117117     endif
    118 
     118     
    119119     dTdp  = 1/(rho_n*cp_n+rho_v*cp_v)
    120120     dPndp = 1/(1d0+m_n/m_v*rho_v/rho_n)

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

    r1525 r1714  
    33  use infotrac, only: nqtot
    44  use radinc_h,      only: NAERKIND
    5   use radcommon_h,   only: L_NSPECTI, L_NSPECTV, sigma
     5  use radcommon_h,   only: L_NSPECTI, L_NSPECTV, sigma, glat
    66  use watercommon_h, only: mH2O
    77  use ioipsl_getincom, only: getin
    … …  
    1111  use inifis_mod, only: inifis
    1212  use comcstfi_mod
     13  use mod_grid_phy_lmdz, only : regular_lonlat
     14  use regular_lonlat_mod, only: init_regular_lonlat
     15  use physics_distribution_mod, only: init_physics_distribution
     16  use regular_lonlat_mod, only: init_regular_lonlat
     17  use mod_interface_dyn_phys, only: init_interface_dyn_phys
     18  use geometry_mod, only: init_geometry
     19  use dimphy, only : init_dimphy
     20
    1321  implicit none
    1422
    … …  
    2735  !     Authors
    2836  !     -------
    29   !     R. Wordsworth
     37  !     - R. Wordsworth
     38  !     - updated by M. Turbet (June 2017)
    3039  !
    3140  !==================================================================
    … …  
    4958  real psurf,psurf_n,tsurf     
    5059
    51   real emis, albedo
     60  real emis, albedo, albedo_equivalent
     61  real albedo_wv(1,L_NSPECTV)
    5262
    5363  real muvar(llm+1)
    … …  
    5868  real fluxtop_lw    ! outgoing LW flux to space (W/m2)
    5969  real fluxsurf_sw   ! incident SW flux to surf (W/m2)
     70  real fluxsurfabs_sw ! absorbed SW flux by the surf (W/m2)
    6071  real fluxabs_sw    ! SW flux absorbed by planet (W/m2)
    6172  real fluxtop_dn    ! incident top of atmosphere SW flux (W/m2)
    6273
    6374  ! not used
    64   real reffrad(llm,naerkind)
    65   real nueffrad(llm,naerkind)
    6675  real cloudfrac(llm)
    6776  real totcloudfrac
    … …  
    7887
    7988  character*20,allocatable :: nametrac(:)   ! name of the tracer (no need for adv trac common)
    80 
     89 
     90  real :: latitude(1), longitude(1), cell_area(1), phisfi(1)
     91 
    8192  ! --------------
    8293  ! Initialisation
    … …  
    8596  pi=2.E+0*asin(1.E+0)
    8697
    87   reffrad(:,:)  = 0.0
    88   nueffrad(:,:) = 0.0
    8998  cloudfrac(:)  = 0.0
    9099  totcloudfrac  = 0.0
    … …  
    97106  ALLOCATE(mu0(1))
    98107  ALLOCATE(fract(1))
    99 
    100 
     108  ALLOCATE(glat(1))
    101109
    102110  !  Load parameters from "run.def"
    … …  
    118126     close(90)
    119127  endif
    120 
     128 
    121129! now, run.def is needed anyway. so we create a dummy temporary one
    122130! for ioipsl to work. if a run.def is already here, stop the
    … …  
    136144     call system("echo 'INCLUDEDEF=kcm1d.def' >> run.def")
    137145  endif
    138 
    139   global1d = .false. ! default value
    140   call getin("global1d",global1d)
    141   if(.not.global1d)then
    142      print*,'Error, kcm1d must have global1d=.true. in kcm1d.def!'
    143      stop
    144   end if
    145 
    146   psurf_n=100000. ! default value for psurf
    147   print*,'Dry surface pressure (Pa)?'
    148   call getin("psurf",psurf_n)
    149   write(*,*) " psurf = ",psurf_n
    150 
    151 ! OK. now that run.def has been read once -- any variable is in memory.
    152 ! so we can dump the dummy run.def
    153   call system("rm -rf run.def")
    154 
    155   tsurf=300.0 ! default value for tsurf
    156   print*,'Surface temperature (K)?'
    157   call getin("tref",tsurf)
    158   write(*,*) " tsurf = ",tsurf
    159 
    160   g=10.0 ! default value for g
    161   print*,'Gravity ?'
    162   call getin("g",g)
    163   write(*,*) " g = ",g
    164 
    165   periastr = 1.0
    166   apoastr  = 1.0
    167   print*,'Periastron (AU)?'
    168   call getin("periastr",periastr)
    169   write(*,*) "periastron = ",periastr
    170   dist_star = periastr
    171   fract     = 0.5
    172   print*,'Apoastron (AU)?'
    173   call getin("apoastr",apoastr)
    174   write(*,*) "apoastron = ",apoastr
    175 
    176   albedo=0.2 ! default value for albedo
    177   print*,'Albedo of bare ground?'
    178   call getin("albedo",albedo)
    179   write(*,*) " albedo = ",albedo
    180 
    181   emis=1.0 ! default value for emissivity
    182   PRINT *,'Emissivity of bare ground ?'
    183   call getin("emis",emis)
    184   write(*,*) " emis = ",emis
    185 
    186   pceil=100.0 ! Pascals
    187   PRINT *,'Ceiling pressure (Pa) ?'
    188   call getin("pceil",pceil)
    189   write(*,*) " pceil = ", pceil
    190 
    191   mugaz=0.
    192   cpp=0.
    193 
    194   check_cpp_match = .false.
    195   call getin("check_cpp_match",check_cpp_match)
    196   if (check_cpp_match) then
    197      print*,"In 1D modeling, check_cpp_match is supposed to be F"
    198      print*,"Please correct callphys.def"
    199      stop
    200   endif
    201 
    202   call su_gases
    203   call calc_cpp_mugaz
    204 
    205   call inifis(1,llm,0,1,86400.0,1,1.0,&
    206               (/0.0/),(/0.0/),(/1.0/),&
    207               rad,g,r,cpp)
     146   
     147! check if we are going to run with or without tracers
     148  write(*,*) "Run with or without tracer transport ?"
     149  tracer=.false. ! default value
     150  call getin("tracer",tracer)
     151  write(*,*) " tracer = ",tracer
     152
     153
    208154
    209155  ! Tracer initialisation
    … …  
    232178           ! minimal version, just read in the tracer names, 1 per line
    233179           read(90,*,iostat=ierr) nametrac(iq)
     180           write(*,*) 'tracer here is', nametrac(iq)
    234181           if (ierr.ne.0) then
    235182              write(*,*) 'kcm1d: error reading tracer names...'
    … …  
    242189
    243190  endif
    244 
     191 
     192
     193
     194
     195  global1d = .false. ! default value
     196  call getin("global1d",global1d)
     197  write(*,*) " global1d = ",global1d
     198  if(.not.global1d)then
     199     print*,'Error, kcm1d must have global1d=.true. in kcm1d.def!'
     200     stop
     201  end if
     202
     203  psurf_n=100000. ! default value for psurf
     204  print*,'Dry surface pressure (Pa)?'
     205  call getin("psurf",psurf_n)
     206  write(*,*) " psurf = ",psurf_n
     207
     208! OK. now that run.def has been read once -- any variable is in memory.
     209! so we can dump the dummy run.def
     210  call system("rm -rf run.def")
     211
     212  tsurf=300.0 ! default value for tsurf
     213  print*,'Surface temperature (K)?'
     214  call getin("tref",tsurf)
     215  write(*,*) " tsurf = ",tsurf
     216
     217  g=10.0 ! default value for g
     218  print*,'Gravity ?'
     219  call getin("g",g)
     220  write(*,*) " g = ",g
     221  glat(1)=g
     222
     223  periastr = 1.0
     224  apoastr  = 1.0
     225  print*,'Periastron (AU)?'
     226  call getin("periastr",periastr)
     227  write(*,*) "periastron = ",periastr
     228  dist_star = periastr
     229  fract     = 0.5
     230  print*,'Apoastron (AU)?'
     231  call getin("apoastr",apoastr)
     232  write(*,*) "apoastron = ",apoastr
     233
     234  albedo=0.2 ! default value for albedo
     235  print*,'Albedo of bare ground?'
     236  call getin("albedo",albedo)
     237  write(*,*) " albedo = ",albedo
     238  do iw=1,L_NSPECTV
     239     albedo_wv(1,iw)=albedo
     240  enddo
     241
     242  emis=1.0 ! default value for emissivity
     243  PRINT *,'Emissivity of bare ground ?'
     244  call getin("emis",emis)
     245  write(*,*) " emis = ",emis
     246
     247  pceil=100.0 ! Pascals
     248  PRINT *,'Ceiling pressure (Pa) ?'
     249  call getin("pceil",pceil)
     250  write(*,*) " pceil = ", pceil
     251
     252  check_cpp_match = .false.
     253  call getin("check_cpp_match",check_cpp_match)
     254  if (check_cpp_match) then
     255     print*,"In 1D modeling, check_cpp_match is supposed to be F"
     256     print*,"Please correct callphys.def"
     257     stop
     258  endif
     259
     260
     261!!! GEOGRAPHICAL INITIALIZATIONS
     262     !!! AREA. useless in 1D
     263  cell_area(1)=1.E+0 !JL+EM to have access to the area in the diagfi.nc files.
     264     !!! GEOPOTENTIAL. useless in 1D because control by surface pressure
     265  phisfi(1)=0.E+0
     266     !!! LATITUDE. can be set.
     267  latitude=0 ! default value for latitude
     268  PRINT *,'latitude (in degrees) ?'
     269  call getin("latitude",latitude)
     270  write(*,*) " latitude = ",latitude
     271  latitude=latitude*pi/180.E+0
     272     !!! LONGITUDE. useless in 1D.
     273  longitude=0.E+0
     274  longitude=longitude*pi/180.E+0
     275
     276  rad=6400000.
     277  !rad = -99999.
     278  !PRINT *,'PLANETARY RADIUS in m ?'
     279  !call getin("rad",rad)
     280  ! Planetary  radius is needed to compute shadow of the rings
     281  !IF (rad.eq.-99999. .and. rings_shadow .eqv. .true.) THEN
     282  !   PRINT *,"STOP. I NEED rad IN RCM1D.DEF."
     283  !   STOP
     284  !ELSE
     285  !   PRINT *,"--> rad = ",rad
     286  !ENDIF
     287
     288
     289
     290  call init_physics_distribution(regular_lonlat,4,1,1,1,nlayer,1)
     291  call init_interface_dyn_phys
     292  CALL init_regular_lonlat(1,1,longitude,latitude,(/0.,0./),(/0.,0./))
     293  call init_geometry(1,longitude,latitude,(/0.,0.,0.,0./),(/0.,0.,0.,0./),cell_area)
     294  call init_dimphy(1,nlayer) ! Initialize dimphy module
     295
     296!!! CALL INIFIS
     297!!! - read callphys.def
     298!!! - calculate sine and cosine of longitude and latitude
     299!!! - calculate mugaz and cp
     300!!! NB: some operations are being done dummily in inifis in 1D
     301!!! - physical constants: nevermind, things are done allright below
     302!!! - physical frequency: nevermind, in inifis this is a simple print
     303
     304  cpp=1.d-7 !JL because we divide by cpp in inifis, there may be a more elegant solution
     305  CALL inifis(1,llm,nq,0,86400.0,1,1.0,latitude,longitude,cell_area,rad,g,r,cpp)
     306
     307  !write(*,*) 'BASE 1'
     308  !write(*,*) 1,llm,nq,0,86400.0,1,1.0,latitude,longitude,cell_area,rad,g,r,cpp
    245309
    246310  do iq=1,nq
    … …  
    249313     enddo
    250314  enddo
    251 
     315 
    252316  do iq=1,nq
    253317     qsurf(iq) = 0.
    … …  
    265329  ! ---------
    266330  !do
    267      psurf = psurf_n
     331  psurf = psurf_n
    268332
    269333     !    Create vertical profiles
    270      call kcmprof_fn(nlayer,psurf,qsurf(1),tsurf,    &
    271           tstrat,play,plev,zlay,temp,q(:,1),muvar(1))
     334  call kcmprof_fn(nlayer,psurf,qsurf(1),tsurf,    &
     335       tstrat,play,plev,zlay,temp,q(:,1),muvar(1))
     336         
     337
     338  !write(*,*) 'BASE 2'
     339  !write(*,*) nlayer,psurf,qsurf(1),tsurf,    &
     340  !           tstrat
     341  !write(*,*) play
     342  !write(*,*) plev
     343  !write(*,*) zlay
     344  !write(*,*) temp
     345  !write(*,*) q(:,1),muvar(1)
     346
     347         
     348               
     349     call iniaerosol()
     350     
     351
     352     
    272353
    273354     !    Run radiative transfer
    274      call callcorrk(1,nlayer,q,nq,qsurf,      &
    275           albedo,emis,mu0,plev,play,temp,                    &
    276           tsurf,fract,dist_star,aerosol,muvar,         &
    277           dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
    278           fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,reffrad,nueffrad,tau_col,  &
     355     call callcorrk(1,nlayer,q,nq,qsurf,                  &
     356          albedo_wv,albedo_equivalent,                    &
     357          emis,mu0,plev,play,temp,                        &
     358          tsurf,fract,dist_star,aerosol,muvar,            &
     359          dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,              &
     360          fluxsurfabs_sw,fluxtop_lw,                      &
     361          fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,tau_col,    &
    279362          cloudfrac,totcloudfrac,.false.,firstcall,lastcall)
     363
     364     !write(*,*) 'BASE 3'
     365     !write(*,*) 1,nlayer,'A',q,'B',nq,'C',qsurf,                      &
     366     !     albedo_wv,'D',albedo_equivalent,'E',                    &
     367     !     emis,'F',mu0,'G',plev,'H',play,'I',temp,'J',                        &
     368     !     tsurf,'K',fract,'L',dist_star,'M',aerosol,'N',muvar,'O',            &
     369     !     dtlw,'P',dtsw,'Q',fluxsurf_lw,'R',fluxsurf_sw,'S',              &
     370     !     fluxsurfabs_sw,'T',fluxtop_lw,'U',                      &
     371     !     fluxabs_sw,'V',fluxtop_dn,'W',OLR_nu,'X',OSR_nu,'Y',tau_col,'Z',    &
     372     !     cloudfrac,'A1',totcloudfrac,'A2',.false.,'A3',firstcall,'A4',lastcall
     373
     374
     375
     376
     377
     378
    280379
    281380     !    Iterate stratospheric temperature
    … …  
    302401  !end do
    303402
     403
    304404  ! Run radiative transfer one last time to get OLR,OSR
    305405  firstcall=.false.
    306406  lastcall=.true.
    307   call callcorrk(1,nlayer,q,nq,qsurf,      &
    308        albedo,emis,mu0,plev,play,temp,                    &
    309        tsurf,fract,dist_star,aerosol,muvar,         &
    310        dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
    311        fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,             &
    312        reffrad,nueffrad,tau_col,  &
     407  call callcorrk(1,nlayer,q,nq,qsurf,                          &
     408       albedo_wv,albedo_equivalent,emis,mu0,plev,play,temp,    &
     409       tsurf,fract,dist_star,aerosol,muvar,                    &
     410       dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxsurfabs_sw,       &
     411       fluxtop_lw, fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,        &
     412       tau_col,                                                &
    313413       cloudfrac,totcloudfrac,.false.,firstcall,lastcall)
     414
     415
     416     !write(*,*) 'BASE 4'
     417     !write(*,*) 1,nlayer,'A',q,'B',nq,'C',qsurf,                      &
     418     !     albedo_wv,'D',albedo_equivalent,'E',                    &
     419     !     emis,'F',mu0,'G',plev,'H',play,'I',temp,'J',                        &
     420     !     tsurf,'K',fract,'L',dist_star,'M',aerosol,'N',muvar,'O',            &
     421     !     dtlw,'P',dtsw,'Q',fluxsurf_lw,'R',fluxsurf_sw,'S',              &
     422     !     fluxsurfabs_sw,'T',fluxtop_lw,'U',                      &
     423     !     fluxabs_sw,'V',fluxtop_dn,'W',OLR_nu,'X',OSR_nu,'Y',tau_col,'Z',    &
     424     !     cloudfrac,'A1',totcloudfrac,'A2',.false.,'A3',firstcall,'A4',lastcall
     425
     426
     427
     428
     429
    314430
    315431

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

    r1403 r1714  
    183183  ! log pressure grid spacing (constant)
    184184  dlogp = -(log(Pn(1)+Pv(1))-log(ptop))/(nlay-1)
    185 
     185 
    186186  call gradients_kcm(profil_flag(1),rho_v(1),rho_n(1),Tsurf,dTdp,dPvdp,dPndp)
    187187  if(verbose)then
Note: See TracChangeset for help on using the changeset viewer.