Changeset 1483 for trunk/LMDZ.GENERIC


Ignore:
Timestamp:
Oct 14, 2015, 8:56:32 PM (9 years ago)
Author:
mturbet
Message:

cleanup of callcorrk

Location:
trunk/LMDZ.GENERIC
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.GENERIC/README

    r1482 r1483  
    10911091- CO2 Ice Albedo North/South dichotomy has been removed. CO2 ice Albedo is no longer stocked in start file. You can now
    10921092prescribe CO2 ice albedo in callphys.def with "albedoco2ice". Its default value is 0.5.
     1093
     1094== 14/10/2015 == MT
     1095- Cleanup of callcorrk.F90
  • trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90

    r1482 r1483  
    1313      use watercommon_h
    1414      use datafile_mod, only: datadir
    15 !      use ioipsl_getincom
    1615      use ioipsl_getincom_p
    1716      use gases_h
     
    2120      use comcstfi_mod, only: pi, mugaz, cpp
    2221      use callkeys_mod, only: varactive,diurnal,tracer,water,nosurf,varfixed,satval,    &
    23                 kastprof,strictboundcorrk,specOLR,CLFvarying
     22                              kastprof,strictboundcorrk,specOLR,CLFvarying
    2423
    2524      implicit none
     
    4443!     Layer #1 is the layer near the ground.
    4544!     Layer #nlayer is the layer at the top.
    46 
    47       INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
    48       INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
    49       REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
    50       integer,intent(in) :: nq ! number of tracers
    51       REAL,INTENT(IN) :: qsurf(ngrid,nq) ! tracer on surface (kg.m-2)
    52       REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV)   ! Spectral SW albedo.MT2015.
    53       REAL,INTENT(IN) :: emis(ngrid)     ! LW emissivity
    54       real,intent(in) :: mu0(ngrid) ! cosine of sun incident angle
    55       REAL,INTENT(IN) :: pplev(ngrid,nlayer+1)  ! inter-layer pressure (Pa)
    56       REAL,INTENT(IN) :: pplay(ngrid,nlayer)    ! mid-layer pressure (Pa)
    57       REAL,INTENT(IN) :: pt(ngrid,nlayer)  ! air temperature (K)
    58       REAL,INTENT(IN) :: tsurf(ngrid)        ! surface temperature (K)
    59       REAL,INTENT(IN) :: fract(ngrid)        ! fraction of day
    60       REAL,INTENT(IN) :: dist_star           ! distance star-planet (AU)
    61       REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol tau (kg/kg)
    62       real,intent(in) :: muvar(ngrid,nlayer+1)
    63       REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! heating rate (K/s) due to LW
    64       REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! heating rate (K/s) due to SW
    65       REAL,INTENT(OUT) :: fluxsurf_lw(ngrid)   ! incident LW flux to surf (W/m2)
    66       REAL,INTENT(OUT) :: fluxsurf_sw(ngrid)   ! incident SW flux to surf (W/m2)
    67       REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid)   ! absorbed SW flux by the surface (W/m2). By MT2015.
    68       REAL,INTENT(OUT) :: fluxtop_lw(ngrid)    ! outgoing LW flux to space (W/m2)
    69       REAL,INTENT(OUT) :: fluxabs_sw(ngrid)    ! SW flux absorbed by planet (W/m2)
    70       REAL,INTENT(OUT) :: fluxtop_dn(ngrid)    ! incident top of atmosphere SW flux (W/m2)
    71       REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
    72       REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
    73       REAL,INTENT(OUT) :: tau_col(ngrid) ! diagnostic from aeropacity
    74       REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
    75 !     for H2O cloud fraction in aeropacity
    76       real,intent(in) :: cloudfrac(ngrid,nlayer)
    77       real,intent(out) :: totcloudfrac(ngrid)
     45!-----------------------------------------------------------------------
     46
     47
     48      ! INPUT
     49      INTEGER,INTENT(IN) :: ngrid                  ! Number of atmospheric columns.
     50      INTEGER,INTENT(IN) :: nlayer                 ! Number of atmospheric layers.
     51      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)       ! Tracers (kg/kg_of_air).
     52      INTEGER,INTENT(IN) :: nq                     ! Number of tracers.
     53      REAL,INTENT(IN) :: qsurf(ngrid,nq)           ! Tracers on surface (kg.m-2).
     54      REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV)   ! Spectral Short Wavelengths Albedo. By MT2015
     55      REAL,INTENT(IN) :: emis(ngrid)               ! Long Wave emissivity.
     56      REAL,INTENT(IN) :: mu0(ngrid)                ! Cosine of sun incident angle.
     57      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1)     ! Inter-layer pressure (Pa).
     58      REAL,INTENT(IN) :: pplay(ngrid,nlayer)       ! Mid-layer pressure (Pa).
     59      REAL,INTENT(IN) :: pt(ngrid,nlayer)          ! Air temperature (K).
     60      REAL,INTENT(IN) :: tsurf(ngrid)              ! Surface temperature (K).
     61      REAL,INTENT(IN) :: fract(ngrid)              ! Fraction of day.
     62      REAL,INTENT(IN) :: dist_star                 ! Distance star-planet (AU).
     63      REAL,INTENT(IN) :: muvar(ngrid,nlayer+1)
     64      REAL,INTENT(IN) :: cloudfrac(ngrid,nlayer)   ! Fraction of clouds (%).
    7865      logical,intent(in) :: clearsky
    79       logical,intent(in) :: firstcall ! signals first call to physics
    80       logical,intent(in) :: lastcall ! signals last call to physics
    81 
    82 !     Globally varying aerosol optical properties on GCM grid
    83 !     Not needed everywhere so not in radcommon_h
     66      logical,intent(in) :: firstcall              ! Signals first call to physics.
     67      logical,intent(in) :: lastcall               ! Signals last call to physics.
     68     
     69      ! OUTPUT
     70      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau (kg/kg).
     71      REAL,INTENT(OUT) :: dtlw(ngrid,nlayer)             ! Heating rate (K/s) due to LW radiation.
     72      REAL,INTENT(OUT) :: dtsw(ngrid,nlayer)             ! Heating rate (K/s) due to SW radiation.
     73      REAL,INTENT(OUT) :: fluxsurf_lw(ngrid)             ! Incident LW flux to surf (W/m2).
     74      REAL,INTENT(OUT) :: fluxsurf_sw(ngrid)             ! Incident SW flux to surf (W/m2)
     75      REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid)          ! Absorbed SW flux by the surface (W/m2). By MT2015.
     76      REAL,INTENT(OUT) :: fluxtop_lw(ngrid)              ! Outgoing LW flux to space (W/m2).
     77      REAL,INTENT(OUT) :: fluxabs_sw(ngrid)              ! SW flux absorbed by the planet (W/m2).
     78      REAL,INTENT(OUT) :: fluxtop_dn(ngrid)              ! Incident top of atmosphere SW flux (W/m2).
     79      REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI)        ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1).
     80      REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV)        ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1).
     81      REAL,INTENT(OUT) :: tau_col(ngrid)                 ! Diagnostic from aeropacity.
     82      REAL,INTENT(OUT) :: albedo_equivalent(ngrid)       ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
     83      REAL,INTENT(OUT) :: totcloudfrac(ngrid)            ! Column Fraction of clouds (%).
     84     
     85     
     86     
     87     
     88
     89      ! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h.   
    8490      REAL :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind)
    8591      REAL :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
    8692      REAL :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
    87 
    8893      REAL :: QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind)
    8994      REAL :: omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
     
    9398!      REAL :: omegaREFir3d(ngrid,nlayer,naerkind) ! not sure of the point of these...
    9499
    95       REAL,ALLOCATABLE,SAVE :: reffrad(:,:,:) ! aerosol effective radius (m)
     100      REAL,ALLOCATABLE,SAVE :: reffrad(:,:,:)  ! aerosol effective radius (m)
    96101      REAL,ALLOCATABLE,SAVE :: nueffrad(:,:,:) ! aerosol effective variance
    97102!$OMP THREADPRIVATE(reffrad,nueffrad)
     
    99104!-----------------------------------------------------------------------
    100105!     Declaration of the variables required by correlated-k subroutines
    101 !     Numbered from top to bottom unlike in the GCM!
     106!     Numbered from top to bottom (unlike in the GCM)
     107!-----------------------------------------------------------------------
    102108
    103109      REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
    104110      REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
    105111
    106 !    Optical values for the optci/cv subroutines
     112      ! Optical values for the optci/cv subroutines
    107113      REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
    108114      REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     
    118124      REAL*8 tauaero(L_LEVELS+1,naerkind)
    119125      REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
    120       real*8 nfluxoutv_nu(L_NSPECTV) ! outgoing band-resolved VI flux at TOA (W/m2)
    121       real*8 nfluxtopi_nu(L_NSPECTI) ! net band-resolved IR flux at TOA (W/m2)
    122       real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! for 1D diagnostic
     126      REAL*8 nfluxoutv_nu(L_NSPECTV)                 ! Outgoing band-resolved VI flux at TOA (W/m2).
     127      REAL*8 nfluxtopi_nu(L_NSPECTI)                 ! Net band-resolved IR flux at TOA (W/m2).
     128      REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI)         ! For 1D diagnostic.
    123129      REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
    124130      REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
    125131      REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
    126132      REAL*8 albi,acosz
    127       REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo.
    128 
    129       INTEGER ig,l,k,nw,iaer,irad
    130       INTEGER icount
     133      REAL*8 albv(L_NSPECTV)                         ! Spectral Visible Albedo.
     134
     135      INTEGER ig,l,k,nw,iaer
    131136
    132137      real szangle
     
    136141      real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
    137142      real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
    138 
    139       real*8 qvar(L_LEVELS)          ! mixing ratio of variable component (mol/mol)
    140 
    141 !     Local aerosol optical properties for each column on RADIATIVE grid
     143      real*8 qvar(L_LEVELS)                    ! Mixing ratio of variable component (mol/mol).
     144
     145      ! Local aerosol optical properties for each column on RADIATIVE grid.
    142146      real*8,save ::  QXVAER(L_LEVELS+1,L_NSPECTV,naerkind)
    143147      real*8,save ::  QSVAER(L_LEVELS+1,L_NSPECTV,naerkind)
     
    155159
    156160
    157 !     Misc.
    158       logical nantest
    159       real*8  tempv(L_NSPECTV)
    160       real*8  tempi(L_NSPECTI)
     161      ! Miscellaneous :
    161162      real*8  temp,temp1,temp2,pweight
    162163      character(len=10) :: tmp1
    163164      character(len=10) :: tmp2
    164165
    165 !     for fixed water vapour profiles
     166      ! For fixed water vapour profiles.
    166167      integer i_var
    167168      real RH
    168169      real*8 pq_temp(nlayer)
     170! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
    169171      real ptemp, Ttemp, qsat
    170172
    171 !      real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
    172 
    173       !real ptime, pday
    174173      logical OLRz
    175174      real*8 NFLUXGNDV_nu(L_NSPECTV)
    176175
    177 
    178       ! for weird cloud test
    179       real pqtest(ngrid,nlayer,nq)
    180 
    181       real maxrad, minrad
    182            
    183       real,external :: CBRT
    184 
    185 !     included by RW for runaway greenhouse 1D study
     176      ! Included by RW for runaway greenhouse 1D study.
    186177      real vtmp(nlayer)
    187178      REAL*8 muvarrad(L_LEVELS)
    188179     
    189 !     included by MT for albedo calculations.     
     180      ! Included by MT for albedo calculations.     
    190181      REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
    191182
     183
     184
    192185!===============================================================
    193 !     Initialization on first call
     186!           I.a Initialization on first call
     187!===============================================================
     188
    194189
    195190      qxvaer(:,:,:)=0.0
     
    203198      if(firstcall) then
    204199
    205          !!! ALLOCATED instances are necessary because of CLFvarying
    206          !!! strategy to call callcorrk twice in physiq...
     200         !!! ALLOCATED instances are necessary because of CLFvarying (strategy to call callcorrk twice in physiq...)
    207201         IF(.not.ALLOCATED(QREFvis3d)) ALLOCATE(QREFvis3d(ngrid,nlayer,naerkind))
    208202         IF(.not.ALLOCATED(QREFir3d)) ALLOCATE(QREFir3d(ngrid,nlayer,naerkind))
     
    219213         call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
    220214         
     215         
    221216!--------------------------------------------------
    222 !     set up correlated k
     217!             Set up correlated k
     218!--------------------------------------------------
     219
     220
    223221         print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
    224222         call getin_p("corrkdir",corrkdir)
     
    229227         banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
    230228
    231          call setspi            ! basic infrared properties
    232          call setspv            ! basic visible properties
    233          call sugas_corrk       ! set up gaseous absorption properties
    234          call suaer_corrk       ! set up aerosol optical properties
     229         call setspi            ! Basic infrared properties.
     230         call setspv            ! Basic visible properties.
     231         call sugas_corrk       ! Set up gaseous absorption properties.
     232         call suaer_corrk       ! Set up aerosol optical properties.
    235233
    236234
     
    244242
    245243         if (ngrid.eq.1) then
    246            PRINT*, 'Simulate global averaged conditions ?'
    247            global1d = .false. ! default value
    248            call getin_p("global1d",global1d)
    249            write(*,*) "global1d = ",global1d
    250            ! Test of incompatibility:
    251            ! if global1d is true, there should not be any diurnal cycle
    252            if (global1d.and.diurnal) then
    253             print*,'if global1d is true, diurnal must be set to false'
    254             stop
    255            endif
    256 
    257            if (global1d) then
    258              PRINT *,'Solar Zenith angle (deg.) ?'
    259              PRINT *,'(assumed for averaged solar flux S/4)'
    260              szangle=60.0  ! default value
    261              call getin_p("szangle",szangle)
    262              write(*,*) "szangle = ",szangle
    263            endif
     244            PRINT*, 'Simulate global averaged conditions ?'
     245            global1d = .false. ! default value
     246            call getin_p("global1d",global1d)
     247            write(*,*) "global1d = ",global1d
     248            
     249            ! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
     250            if (global1d.and.diurnal) then
     251               print*,'if global1d is true, diurnal must be set to false'
     252               stop
     253            endif
     254
     255            if (global1d) then
     256               PRINT *,'Solar Zenith angle (deg.) ?'
     257               PRINT *,'(assumed for averaged solar flux S/4)'
     258               szangle=60.0  ! default value
     259               call getin_p("szangle",szangle)
     260               write(*,*) "szangle = ",szangle
     261            endif
    264262         endif
    265263
    266264      end if ! of if (firstcall)
    267265
     266
    268267!=======================================================================
    269 !     Initialization on every call   
    270 
     268!          I.b  Initialization on every call   
     269!=======================================================================
     270 
     271 
    271272!--------------------------------------------------
    272273!     Effective radius and variance of the aerosols
     274!--------------------------------------------------
     275
    273276      do iaer=1,naerkind
    274277
    275          if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! treat condensed co2 particles.
     278         if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! Treat condensed co2 particles.
    276279            call co2_reffrad(ngrid,nlayer,nq,pq,reffrad(1,1,iaero_co2))
    277280            print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
    278281            print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
    279282         end if
    280          if ((iaer.eq.iaero_h2o).and.water) then ! treat condensed water particles. to be generalized for other aerosols
     283         
     284         if ((iaer.eq.iaero_h2o).and.water) then ! Treat condensed water particles. To be generalized for other aerosols ...
    281285            call h2o_reffrad(ngrid,nlayer,pq(1,1,igcm_h2o_ice),pt, &
    282286                             reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
     
    284288            print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
    285289         endif
     290         
    286291         if(iaer.eq.iaero_dust)then
    287292            call dust_reffrad(ngrid,nlayer,reffrad(1,1,iaero_dust))
    288293            print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
    289294         endif
     295         
    290296         if(iaer.eq.iaero_h2so4)then
    291297            call h2so4_reffrad(ngrid,nlayer,reffrad(1,1,iaero_h2so4))
    292298            print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
    293299         endif
     300         
    294301          if(iaer.eq.iaero_back2lay)then
    295302            call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev)
    296303         endif
    297      end do !iaer=1,naerkind
    298 
    299 
    300 !     how much light we get
     304         
     305     end do !iaer=1,naerkind.
     306
     307
     308      ! How much light do we get ?
    301309      do nw=1,L_NSPECTV
    302310         stel(nw)=stellarf(nw)/(dist_star**2)
    303311      end do
    304312
     313      ! Get 3D aerosol optical properties.
    305314      call aeroptproperties(ngrid,nlayer,reffrad,nueffrad,         &
    306315           QVISsQREF3d,omegaVIS3d,gVIS3d,                          &
    307316           QIRsQREF3d,omegaIR3d,gIR3d,                             &
    308            QREFvis3d,QREFir3d)                                     ! get 3D aerosol optical properties
    309 
     317           QREFvis3d,QREFir3d)                                     
     318
     319      ! Get aerosol optical depths.
    310320      call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol,      &
    311321           reffrad,QREFvis3d,QREFir3d,                             &
    312            tau_col,cloudfrac,totcloudfrac,clearsky)                ! get aerosol optical depths
     322           tau_col,cloudfrac,totcloudfrac,clearsky)               
    313323         
    314 !-----------------------------------------------------------------------
    315 !     Starting Big Loop over every GCM column
    316       do ig=1,ngrid
     324
     325
     326!-----------------------------------------------------------------------   
     327      do ig=1,ngrid ! Starting Big Loop over every GCM column
     328!-----------------------------------------------------------------------
     329
    317330
    318331!=======================================================================
    319 !     Transformation of the GCM variables
    320 
    321 !-----------------------------------------------------------------------
    322 !     Aerosol optical properties Qext, Qscat and g
    323 !     The transformation in the vertical is the same as for temperature
     332!              II.  Transformation of the GCM variables
     333!=======================================================================
     334
     335
     336!-----------------------------------------------------------------------
     337!    Aerosol optical properties Qext, Qscat and g.
     338!    The transformation in the vertical is the same as for temperature.
     339!-----------------------------------------------------------------------
    324340           
    325 !     shortwave
     341           
    326342            do iaer=1,naerkind
    327                DO nw=1,L_NSPECTV
     343               ! Shortwave.
     344               do nw=1,L_NSPECTV
     345               
    328346                  do l=1,nlayer
    329347
     
    349367                     gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
    350368
    351                   end do
     369                  end do ! nlayer
    352370
    353371                  qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
     
    360378                  gvaer(2*nlayer+1,nw,iaer)=0.
    361379
    362                end do
    363 
    364 !     longwave
    365                DO nw=1,L_NSPECTI
     380               end do ! L_NSPECTV
     381             
     382               do nw=1,L_NSPECTI
     383                  ! Longwave
    366384                  do l=1,nlayer
    367385
     
    387405                     giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
    388406
    389                   end do
     407                  end do ! nlayer
    390408
    391409                  qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
     
    398416                  giaer(2*nlayer+1,nw,iaer)=0.
    399417
    400                end do
    401             end do
    402 
    403             ! test / correct for freaky s. s. albedo values
     418               end do ! L_NSPECTI
     419               
     420            end do ! naerkind
     421
     422            ! Test / Correct for freaky s. s. albedo values.
    404423            do iaer=1,naerkind
    405424               do k=1,L_LEVELS+1
     
    427446                  end do
    428447
    429                end do
    430             end do
     448               end do ! L_LEVELS
     449            end do ! naerkind
    431450
    432451!-----------------------------------------------------------------------
    433452!     Aerosol optical depths
     453!-----------------------------------------------------------------------
    434454           
    435455         do iaer=1,naerkind     ! a bug was here           
     
    437457               
    438458               pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/   &
    439                         (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
    440 
     459                       (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
    441460               temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
    442 
    443461               tauaero(2*k+2,iaer)=max(temp*pweight,0.d0)
    444462               tauaero(2*k+3,iaer)=max(temp-tauaero(2*k+2,iaer),0.d0)
    445 !
     463
    446464            end do
    447465            ! boundary conditions
     
    450468            !tauaero(1,iaer)          = 0.
    451469            !tauaero(L_LEVELS+1,iaer) = 0.
    452          end do
    453 
    454 !     Albedo and emissivity
    455          albi=1-emis(ig)   ! Longwave.
    456          DO nw=1,L_NSPECTV ! Shortwave loop.
     470           
     471         end do ! naerkind
     472
     473         ! Albedo and Emissivity.
     474         albi=1-emis(ig)   ! Long Wave.
     475         DO nw=1,L_NSPECTV ! Short Wave loop.
    457476            albv(nw)=albedo(ig,nw)
    458477         ENDDO
    459478
    460          if (nosurf) then
     479         if (nosurf) then ! Case with no surface.
    461480            DO nw=1,L_NSPECTV
    462481               if(albv(nw).gt.0.0) then
     
    468487         endif
    469488
    470       if ((ngrid.eq.1).and.(global1d)) then       ! fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight
     489      if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
    471490         acosz = cos(pi*szangle/180.0)
    472491         print*,'acosz=',acosz,', szangle=',szangle
    473492      else
    474          acosz=mu0(ig)          ! cosine of sun incident angle : 3D simulations or local 1D simulations using latitude
     493         acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
    475494      endif
    476495
    477 !!! JL13: in the following, I changed some indices in the interpolations so that the model rsults are less dependent on the number of layers
    478 !!!    the older verions are commented with the commetn !JL13index
    479 
    480 
    481 !-----------------------------------------------------------------------
    482 !     Water vapour (to be generalised for other gases eventually)
     496!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     497!!! Note by JL13 : In the following, some indices were changed in the interpolations,
     498!!!                so that the model results are less dependent on the number of layers !
     499!!!
     500!!!           ---  The older versions are commented with the comment !JL13index  ---
     501!!!
     502!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     503
     504
     505!-----------------------------------------------------------------------
     506!     Water vapour (to be generalised for other gases eventually ...)
     507!-----------------------------------------------------------------------
    483508     
    484509      if(varactive)then
     
    495520      elseif(varfixed)then
    496521
    497          do l=1,nlayer        ! here we will assign fixed water vapour profiles globally
     522         do l=1,nlayer ! Here we will assign fixed water vapour profiles globally.
    498523            RH = satval * ((pplay(ig,l)/pplev(ig,1) - 0.02) / 0.98)
    499524            if(RH.lt.0.0) RH=0.0
     
    511536            qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2
    512537         end do
     538         
    513539         qvar(1)=qvar(2)
    514540
     
    527553            qvar(k) = 1.0D-7
    528554         end do
    529       end if
     555      end if ! varactive/varfixed
    530556
    531557      if(.not.kastprof)then
    532       ! IMPORTANT: Now convert from kg/kg to mol/mol
     558         ! IMPORTANT: Now convert from kg/kg to mol/mol.
    533559         do k=1,L_LEVELS
    534560            qvar(k) = qvar(k)/(epsi+qvar(k)*(1.-epsi))
     
    537563
    538564!-----------------------------------------------------------------------
    539 !     kcm mode only
     565!     kcm mode only !
     566!-----------------------------------------------------------------------
     567
    540568      if(kastprof)then
    541569
    542          ! initial values equivalent to mugaz
     570         ! Initial values equivalent to mugaz.
    543571         DO l=1,nlayer
    544572            muvarrad(2*l)   = mugaz
     
    549577
    550578            DO l=1,nlayer
    551                muvarrad(2*l)   = muvar(ig,nlayer+2-l)
     579               muvarrad(2*l)   =  muvar(ig,nlayer+2-l)
    552580               muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + &
    553                                 muvar(ig,max(nlayer+1-l,1)))/2
     581                                  muvar(ig,max(nlayer+1-l,1)))/2
    554582            END DO
    555583     
    556584            muvarrad(1) = muvarrad(2)
    557             muvarrad(2*nlayer+1)=muvar(ig,1)
     585            muvarrad(2*nlayer+1) = muvar(ig,1)
    558586
    559587            print*,'Recalculating qvar with VARIABLE epsi for kastprof'
    560588            print*,'Assumes that the variable gas is H2O!!!'
    561589            print*,'Assumes that there is only one tracer'
     590           
    562591            !i_var=igcm_h2o_vap
    563592            i_var=1
     593           
    564594            if(nq.gt.1)then
    565595               print*,'Need 1 tracer only to run kcm1d.e'
    566596               stop
    567597            endif
     598           
    568599            do l=1,nlayer
    569600               vtmp(l)=pq(ig,l,i_var)/(epsi+pq(ig,l,i_var)*(1.-epsi))
     
    580611            print*,'Warning: reducing qvar in callcorrk.'
    581612            print*,'Temperature profile no longer consistent ', &
    582                             'with saturated H2O. qsat=',satval
     613                   'with saturated H2O. qsat=',satval
     614                   
    583615            do k=1,L_LEVELS
    584616               qvar(k) = qvar(k)*satval
     
    594626         muvarrad(1) = muvarrad(2)
    595627         muvarrad(2*nlayer+1)=muvar(ig,1)         
    596       endif
    597      
    598       ! Keep values inside limits for which we have radiative transfer coefficients
    599       if(L_REFVAR.gt.1)then ! there was a bug here!
     628      endif ! if kastprof
     629     
     630      ! Keep values inside limits for which we have radiative transfer coefficients !!!
     631      if(L_REFVAR.gt.1)then ! (there was a bug here)
    600632         do k=1,L_LEVELS
    601633            if(qvar(k).lt.wrefvar(1))then
     
    609641!-----------------------------------------------------------------------
    610642!     Pressure and temperature
     643!-----------------------------------------------------------------------
    611644
    612645      DO l=1,nlayer
     
    618651     
    619652      plevrad(1) = 0.
    620       plevrad(2) = 0.   !! Trick to have correct calculations of fluxes in gflux(i/v).F,
    621                         !! but the pmid levels are not impacted by this change.
     653      plevrad(2) = 0.   !! Trick to have correct calculations of fluxes in gflux(i/v).F, but the pmid levels are not impacted by this change.
    622654
    623655      tlevrad(1) = tlevrad(2)
     
    653685!      tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
    654686
    655       ! test for out-of-bounds pressure
     687      ! Test for out-of-bounds pressure.
    656688      if(plevrad(3).lt.pgasmin)then
    657689         print*,'Minimum pressure is outside the radiative'
     
    664696      endif
    665697
    666       ! test for out-of-bounds temperature
     698      ! Test for out-of-bounds temperature.
    667699      do k=1,L_LEVELS
    668700         if(tlevrad(k).lt.tgasmin)then
     
    733765
    734766!=======================================================================
    735 !     Calling the main radiative transfer subroutines
    736 
    737 
    738          Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed
     767!          III. Calling the main radiative transfer subroutines
     768!=======================================================================
     769
     770
     771         Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
    739772         glat_ig=glat(ig)
    740773
    741774!-----------------------------------------------------------------------
    742 !     Shortwave
    743 
    744          if(fract(ig) .ge. 1.0e-4) then ! only during daylight!
     775!        Short Wave Part
     776!-----------------------------------------------------------------------
     777
     778         if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
    745779            if((ngrid.eq.1).and.(global1d))then
    746780               do nw=1,L_NSPECTV
    747                   stel_fract(nw)= stel(nw)* 0.25 / acosz
    748                                 ! globally averaged = divide by 4
    749                                 ! but we correct for solar zenith angle
     781                  stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle
    750782               end do
    751783            else
     
    753785                  stel_fract(nw)= stel(nw) * fract(ig)
    754786               end do
    755             endif
     787            endif
     788           
    756789            call optcv(dtauv,tauv,taucumv,plevrad,                 &
    757790                 qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero,   &
     
    763796                 fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
    764797
    765          else                          ! during the night, fluxes = 0
     798         else ! During the night, fluxes = 0.
    766799            nfluxtopv       = 0.0d0
    767800            fluxtopvdn      = 0.0d0
     
    789822
    790823!-----------------------------------------------------------------------
    791 !     Longwave
     824!        Long Wave Part
     825!-----------------------------------------------------------------------
    792826
    793827         call optci(plevrad,tlevrad,dtaui,taucumi,                  &
     
    850884             *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(1)))
    851885
    852 !     ---------------------------------------------------------------
    853       end do                    ! end of big loop over every GCM column (ig = 1:ngrid)
     886
     887!-----------------------------------------------------------------------   
     888      end do ! End of big loop over every GCM column.
     889!-----------------------------------------------------------------------
     890
    854891
    855892
    856893!-----------------------------------------------------------------------
    857894!     Additional diagnostics
    858 
    859 !     IR spectral output, for exoplanet observational comparison
    860 
    861 
     895!-----------------------------------------------------------------------
     896
     897      ! IR spectral output, for exoplanet observational comparison
    862898      if(lastcall.and.(ngrid.eq.1))then  ! could disable the 1D output, they are in the diagfi and diagspec... JL12
    863899
    864            print*,'Saving scalar quantities in surf_vals.out...'
    865            print*,'psurf = ', pplev(1,1),' Pa'
    866            open(116,file='surf_vals.out')
    867            write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1),         &
    868                 real(-nfluxtopv),real(nfluxtopi)
    869            close(116)
    870 
    871 !          I am useful, please don`t remove me!
     900         print*,'Saving scalar quantities in surf_vals.out...'
     901         print*,'psurf = ', pplev(1,1),' Pa'
     902         open(116,file='surf_vals.out')
     903         write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1),         &
     904                      real(-nfluxtopv),real(nfluxtopi)
     905         close(116)
     906
     907
     908!          USEFUL COMMENT - Do Not Remove.
     909!
    872910!           if(specOLR)then
    873911!               open(117,file='OLRnu.out')
     
    884922!           endif
    885923
    886 !     OLR vs altitude: do it as a .txt file
    887            OLRz=.false.
    888            if(OLRz)then
    889               print*,'saving IR vertical flux for OLRz...'
    890               open(118,file='OLRz_plevs.out')
    891               open(119,file='OLRz.out')
    892               do l=1,L_NLAYRAD
    893                  write(118,*) plevrad(2*l)
    894                  do nw=1,L_NSPECTI
    895                      write(119,*) fluxupi_nu(l,nw)
    896                   enddo
    897               enddo
    898               close(118)
    899               close(119)
    900            endif
     924           ! OLR vs altitude: do it as a .txt file.
     925         OLRz=.false.
     926         if(OLRz)then
     927            print*,'saving IR vertical flux for OLRz...'
     928            open(118,file='OLRz_plevs.out')
     929            open(119,file='OLRz.out')
     930            do l=1,L_NLAYRAD
     931               write(118,*) plevrad(2*l)
     932               do nw=1,L_NSPECTI
     933                  write(119,*) fluxupi_nu(l,nw)
     934               enddo
     935            enddo
     936            close(118)
     937            close(119)
     938         endif
    901939
    902940      endif
    903941
    904       ! see physiq.F for explanations about CLFvarying. This is temporary.
     942      ! See physiq.F for explanations about CLFvarying. This is temporary.
    905943      if (lastcall .and. .not.CLFvarying) then
    906944        IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
Note: See TracChangeset for help on using the changeset viewer.