Changeset 486 for trunk


Ignore:
Timestamp:
Dec 21, 2011, 10:09:07 AM (13 years ago)
Author:
rwordsworth
Message:

Variable ice timestep added for ice evolution algorithm.
Sedimentation improved to allow consistent dust (L. Kerber).
Bug linked to allocatable matrices removed from rcm1d.F.
Treatment of initial aerosol radii in callcorrk.F90 improved.

Location:
trunk/LMDZ.GENERIC/libf/phystd
Files:
1 deleted
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90

    r305 r486  
    267267               do l=1,nlayer
    268268                  zp=700./pplay(ig,l)
    269                   aerosol(ig,l,1)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) &
     269                  aerosol(ig,l,iaer)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) &
    270270                       *max( exp(.03*(1.-max(zp,1.))) , 1.E-3 ) &
    271271                       *QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)

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

    r471 r486  
    211211!     Effective radius and variance of the aerosols
    212212
    213 !     CO2 ice:
    214          DO l=1,nlayer
    215             DO ig=1,ngrid
    216                reffrad(ig,l,1)  = 1.e-4
    217                nueffrad(ig,l,1) = 0.1
     213         do iaer=1,naerkind
    218214!     these values will change once the microphysics gets to work
    219215!     UNLESS tracer=.false., in which case we should be working with
    220216!     a fixed aerosol layer, and be able to define reffrad in a
    221217!     .def file. To be improved!
    222             ENDDO
    223          ENDDO
    224 
    225 !     H2O ice:
    226          if(naerkind.eq.2)then
    227             DO l=1,nlayer
    228                DO ig=1,ngrid
    229                   reffrad(ig,l,naerkind)  = 1.e-5
    230                   nueffrad(ig,l,naerkind) = 0.1
    231                ENDDO 
    232             ENDDO
    233          endif
     218
     219            if(iaer.eq.1)then ! CO2 ice
     220               do l=1,nlayer
     221                  do ig=1,ngrid
     222                     reffrad(ig,l,iaer)  = 1.e-4
     223                     nueffrad(ig,l,iaer) = 0.1
     224                  enddo
     225               enddo
     226            endif
     227
     228            if(iaer.eq.2)then ! H2O ice
     229               do l=1,nlayer
     230                  do ig=1,ngrid
     231                     reffrad(ig,l,iaer)  = 1.e-5
     232                     nueffrad(ig,l,iaer) = 0.1
     233                  enddo
     234               enddo
     235            endif
     236
     237            if(iaer.eq.3)then ! dust
     238               do l=1,nlayer
     239                  do ig=1,ngrid
     240                     reffrad(ig,l,iaer)  = 1.e-5
     241                     nueffrad(ig,l,iaer) = 0.1
     242                  enddo
     243               enddo
     244            endif
     245 
     246            if(iaer.gt.3)then
     247               print*,'Error in callcorrk, naerkind is too high.'
     248               print*,'The code still needs generalisation to arbitrary'
     249               print*,'aerosol kinds and number.'
     250               call abort
     251            endif
     252
     253         enddo
    234254
    235255         print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
    … …  
    248268         Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed
    249269
    250 
    251270         if((igcm_h2o_vap.eq.0) .and. varactive)then
    252271            print*,'varactive in callcorrk but no h2o_vap tracer.'
    … …  
    268287         enddo
    269288      enddo
    270 
    271289
    272290      if(kastprof)then
    … …  
    426444               end do
    427445
    428 
    429446!     longwave
    430447               DO nw=1,L_NSPECTI
    … …  
    465482               end do
    466483            end do
    467 
    468484
    469485            ! test / correct for freaky s. s. albedo values
    … …  
    515531            !tauaero(L_LEVELS+1,iaer) = 0.
    516532         end do
    517          !print*,'Note changed tauaero BCs in callcorrk!'
    518533
    519534!     Albedo and emissivity
    … …  
    533548         acosz=mu0(ig)          ! cosine of sun incident angle
    534549      endif
    535 
    536550
    537551!-----------------------------------------------------------------------
    … …  
    591605      end do
    592606
    593 
    594 
    595607!-----------------------------------------------------------------------
    596608!     kcm mode only
    597609      if(kastprof)then
    598610
     611         ! initial values equivalent to mugaz
    599612         DO l=1,nlayer
    600613            muvarrad(2*l)   = mugaz
    … …  
    602615         END DO
    603616
    604          do k=1,L_LEVELS
    605             qvar(k) = 0.0
    606          end do
    607          print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!'
     617         !do k=1,L_LEVELS
     618         !   qvar(k) = 0.0
     619         !end do
     620         !print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!'
    608621      endif
    609622
    … …  
    633646         qvar(2*nlayermx+1)=qsurf(ig,i_var)*muvar(ig,1)/mH2O
    634647
    635 
    636       endif
    637 
     648      endif
    638649
    639650      ! Keep values inside limits for which we have radiative transfer coefficients
    … …  
    678689      tmid(L_LEVELS) = tlevrad(L_LEVELS)
    679690
    680 
    681691      ! test for out-of-bounds pressure
    682692      if(plevrad(3).lt.pgasmin)then
    … …  
    695705            print*,'Minimum temperature is outside the radiative'
    696706            print*,'transfer kmatrix bounds, exiting.'
    697             print*,'WARNING, OVERRIDING FOR TEST'
    698             !call abort
     707            !print*,'WARNING, OVERRIDING FOR TEST'
     708            call abort
    699709         elseif(tlevrad(k).gt.tgasmax)then
    700710            print*,'Maximum temperature is outside the radiative'
    701711            print*,'transfer kmatrix bounds, exiting.'
    702             print*,'WARNING, OVERRIDING FOR TEST'
    703             !print*, 'T=',pt
    704             !call abort
     712            !print*,'WARNING, OVERRIDING FOR TEST'
     713            call abort
    705714         endif
    706715      enddo
    … …  
    735744                 tmid,pmid,taugsurf,qvar,muvarrad)
    736745
    737 
    738746            call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
    739747                 acosz,stel_fract,gweight,                         &
    740748                 nfluxtopv,nfluxoutv_nu,nfluxgndv_nu,              &
    741                  !acosz,stel_fract,gweight,nfluxtopv,nfluxgndv_nu,  &
    742749                 fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
    743750
    … …  
    753760         end if
    754761
    755 
    756 
    757 
    758762!-----------------------------------------------------------------------
    759763!     Longwave
    … …  
    846850      if(specOLR)then
    847851        if(ngrid.ne.1)then
    848           !call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu)
    849           !call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W m^-2",3,OSR_nu)
     852          call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu)
     853          call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W m^-2",3,OSR_nu)
    850854        endif
    851855      endif
    … …  
    892896      endif
    893897
    894       !!! see physiq.F for explanations about CLFvarying. This is temporary.
     898      ! see physiq.F for explanations about CLFvarying. This is temporary.
    895899      if (lastcall .and. .not.CLFvarying) then
    896900        IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )

  • trunk/LMDZ.GENERIC/libf/phystd/callkeys.h

    r305 r486  
    3838     &   , hydrology                                                    &
    3939     &   , sourceevol                                                   &
     40     &   , icetstep                                                     &
    4041     &   , albedosnow                                                   &
    4142     &   , maxicethick                                                  &
    … …  
    101102      real tau_relax
    102103      real cloudlvl
     104      real icetstep
     105

  • trunk/LMDZ.GENERIC/libf/phystd/callsedim.F

    r253 r486  
    5353      INTEGER l,ig, iq
    5454
    55       real zqi(ngridmx,nlayermx) ! to locally store tracers
     55      real zqi(ngridmx,nlayermx,nqmx) ! to locally store tracers
    5656      real masse (ngridmx,nlayermx) ! Layer mass (kg.m-2)
    5757      real epaisseur (ngridmx,nlayermx) ! Layer thickness (m)
    … …  
    7979        firstcall=.false.
    8080      ENDIF ! of IF (firstcall)
    81 
     81     
    8282!=======================================================================
    8383!     Preliminary calculation of the layer characteristics
    … …  
    9191      end do
    9292
    93       iq=igcm_h2o_ice
     93!======================================================================
     94! Calculate the transport due to sedimentation for each tracer
     95! [This has been rearranged by L. Kerber to allow the sedimentation
     96!  of general tracers.]
     97 
     98      zqi(1:ngrid,1:nlay,1:nqmx) = 0.
     99      do iq=1,nq
     100       if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_co2_ice)) then   
     101!         (no sedim for gases; co2_ice sedim is done elsewhere)     
    94102
    95 !     The value of q is updated after the other parameterisations
     103! store locally updated tracers
    96104
    97           do l=1,nlay
     105          do l=1,nlay 
    98106            do ig=1,ngrid
    99               ! store locally updated tracers
    100               zqi(ig,l)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
     107              zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
    101108            enddo
    102109          enddo ! of do l=1,nlay
     110       
     111!======================================================================
     112! Sedimentation
     113!======================================================================
     114! Water         
     115          if (water.and.(iq.eq.igcm_h2o_ice)) then
     116             call newsedim(ngrid,nlay,ngrid*nlay,ptimestep,
     117     &            pplev,masse,epaisseur,pt,rfall,rho_q(iq),zqi,wq)
    103118
    104 
    105 
    106 !=======================================================================
    107 !     Calculate the transport due to sedimentation for each tracer
    108 
    109           call newsedim(ngrid,nlay,ngrid*nlay,ptimestep,
    110      &         pplev,masse,epaisseur,pt,rfall,rho_q(iq),zqi,wq)
     119! General Case
     120          else
     121             call newsedim(ngrid,nlay,1,ptimestep,
     122     &            pplev,masse,epaisseur,pt,radius(iq),rho_q(iq),
     123     &            zqi,wq)
     124          endif
    111125
    112126!=======================================================================
    113127!     Calculate the tendencies
     128!======================================================================
    114129
    115130          do ig=1,ngrid
    116131!     Ehouarn: with new way of tracking tracers by name, this is simply
    117             pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
     132            pdqs_sed(ig,iq) = wq(ig,1)/ptimestep
    118133          end do
    119134
    120135          DO l = 1, nlay
    121136            DO ig=1,ngrid
    122               pdqsed(ig,l,iq)=(zqi(ig,l)-
    123      $        (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
     137              pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
     138     &        (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
    124139            ENDDO
    125140          ENDDO
    126 
     141       endif ! of no gases no co2_ice
     142      enddo ! of do iq=1,nq
    127143      return
    128144      end
    129 

  • trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90

    r471 r486  
    187187         enddo
    188188
    189         write(*,*) "condense_co2cloud: i_co2ice=",i_co2ice       
     189        write(*,*) "condense_cloud: i_co2ice=",i_co2ice       
    190190
    191191        if((i_co2ice.lt.1))then
    192            print*,'In condens_co2cloud but no CO2 ice tracer, exiting.'
     192           print*,'In condens_cloud but no CO2 ice tracer, exiting.'
     193           print*,'Still need generalisation to arbitrary species!'
    193194           stop
    194195        endif
    195196
    196197         ccond=cpp/(g*latcond)
    197          print*,'In condens_co2cloud: ccond=',ccond,' latcond=',latcond
     198         print*,'In condens_cloud: ccond=',ccond,' latcond=',latcond
    198199
    199200!          Prepare special treatment if gas is not pure CO2

  • trunk/LMDZ.GENERIC/libf/phystd/inifis.F

    r374 r486  
    472472         write(*,*) " sourceevol = ",sourceevol
    473473
     474         write(*,*) "Ice evolution timestep ?"
     475         icetstep=100.0 ! default value
     476         call getin("icetstep",icetstep)
     477         write(*,*) " icetstep = ",icetstep
     478
    474479         write(*,*) "Snow albedo ?"
    475480         albedosnow=0.5         ! default value

  • trunk/LMDZ.GENERIC/libf/phystd/newsedim.F

    r253 r486  
    7979            PRINT*,'ngrid  =',ngrid
    8080            PRINT*,'ngridmx  =',ngridmx
    81             STOP
     81            STOP 
    8282         ENDIF
    8383         firstcall=.false.
    … …  
    111111c     (stokes law corrected for low pressure by the Cunningham
    112112c     slip-flow correction  according to Rossow (Icarus 36, 1-50, 1978)
    113 
     113       
    114114        do  l=1,nlay
    115115          do ig=1, ngrid
    … …  
    121121            endif 
    122122
    123 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    124 ! TEMPORARY MODIF BY RDW !!!!
    125             !rfall=5.e-6
    126             if(rfall.lt.1.e-7)then
    127                rfall=1.e-7
    128             endif
    129             !if(rfall.gt.5.e-5)then
    130             !   rfall=5.e-5
    131             !endif
    132 
    133 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    134 
    135123            vstokes(ig,l) = b * rho * rfall**2 *
    136124     &      (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rfall)
    … …  
    147135c      va traverser le niveau intercouche l : "dztop" est sa hauteur
    148136c      au dessus de l (m), "ptop" est sa pression (Pa))
    149 
    150137      do  l=1,nlay
    151138        do ig=1, ngrid
    … …  
    154141             Ep=0
    155142             k=0
    156 
     143           w(ig,l) = 0. !! JF+AS ajout initialisation (LK MARS)
    157144c **************************************************************
    158145c            Simple Method
    159              w(ig,l) =
    160      &       (1- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g
    161 cc           write(*,*) 'OK simple method l,w =', l, w(ig,l)
    162 cc           write(*,*) 'OK simple method dztop =', dztop
    163 c **************************************************************
     146cc             w(ig,l) =
     147cc     &       (1- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g
     148cc             write(*,*) 'OK simple method l,w =', l, w(ig,l)
     149cc            write(*,*) 'OK simple method dztop =', dztop
     150           w(ig,l) = 1. - exp(-dztop*g/(r*pt(ig,l)))
     151             !!! Diagnostic: JF. Fix: AS. Date: 05/11
     152             !!! Probleme arrondi avec la quantite ci-dessus
     153             !!! ---> vaut 0 pour -dztop*g/(r*pt(ig,l)) trop petit
     154             !!! ---> dans ce cas on utilise le developpement limite !
     155             !!! ---> exp(-x) = 1 - x lorsque x --> 0 avec une erreur de x^2 / 2 
     156
     157             IF ( w(ig,l) .eq. 0. ) THEN
     158                w(ig,l) = ( dztop*g/(r*pt(ig,l)) ) * pplev(ig,l) / g
     159             ELSE
     160                w(ig,l) = w(ig,l) * pplev(ig,l) / g
     161             ENDIF
     162! LK borrowed simple method from Mars model (AS/JF)
     163
     164!**************************************************************
    164165cccc         Complex method :
    165             if (dztop.gt.epaisseur(ig,l)) then
     166           if (dztop.gt.epaisseur(ig,l)) then
    166167cccc            Cas ou on "epuise" la couche l : On calcule le flux
    167168cccc            Venant de dessus en tenant compte de la variation de Vstokes
    … …  
    176177               enddo
    177178               Ep = Ep - epaisseur(ig,l+k)
    178              end if
    179              ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
    180              w(ig,l) = (pplev(ig,l) -Ptop)/g
     179!             ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
     180             ptop=exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
     181             IF ( ptop .eq. 1. ) THEN
     182                !PRINT*, 'newsedim: exposant trop petit ', ig, l
     183                ptop=pplev(ig,l+k) * ( 1. - (dztop-Ep)*g/(r*pt(ig,l+k)))
     184             ELSE
     185                ptop=pplev(ig,l+k) * ptop
     186             ENDIF
     187
     188             w(ig,l) = (pplev(ig,l) - ptop)/g
     189
     190            endif   !!! complex method
    181191c
    182192cc           write(*,*) 'OK new    method l,w =', l, w(ig,l)
    … …  
    188198cc            if(l.eq.8)write(*,*)'l=8,k, w',k, w(1,l)
    189199c **************************************************************
     200
    190201        end do
    191202      end do
    … …  
    195206c    &                wq(1,6),wq(1,7),pqi(1,6)
    196207
    197 
    198208      RETURN
    199209      END
    200 

  • trunk/LMDZ.GENERIC/libf/phystd/physiq.F90

    r471 r486  
    19171917                  qsurf_hist(ig,igcm_h2o_ice) = &
    19181918                     !qsurf_hist(ig,igcm_h2o_ice) + delta_ice*100.0
    1919                      qsurf_hist(ig,igcm_h2o_ice) + delta_ice*10.0
     1919                     qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep
    19201920
    19211921                  ! if ice has gone -ve, set to zero

  • trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90

    r470 r486  
    66#include "bands.h"
    77
    8 !======================================================================C
     8!======================================================================
    99!
    1010!     RADINC.H    RADiation INCludes
    … …  
    1313!     number of spectral intervals. . .
    1414!
    15 !======================================================================C
     15!======================================================================
    1616
    1717!     RADIATION parameters
    … …  
    6060      integer, parameter :: L_NLEVRAD  = llm+1
    6161
    62       !!!! THIS IS SET IN sugas_corrk
    63       !!!! [use of allocatable arrays] -- AS 12/2011
     62      ! These are set in sugas_corrk
     63      ! [uses allocatable arrays] -- AS 12/2011
    6464      integer :: L_NPREF, L_NTREF, L_REFVAR, L_PINT
    6565
    66 !!! THIS ONE DOES NOT CHANGE MUCH... IT CAN BE REGARDED AS A CONSTANT.
    6766      integer, parameter :: L_NGAUSS  = 17
    6867
    … …  
    7271      integer, parameter :: NAERKIND  = 2
    7372      real,    parameter :: L_TAUMAX  = 35
    74       !integer, parameter :: L_TAUMAX  = 35
    75       !integer, parameter :: L_TAUMAX  = 40
    7673
    7774      ! For Planck function integration:

  • trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F

    r305 r486  
    112112      REAL cloudfrac(ngridmx,nlayermx)
    113113      REAL hice(ngridmx),totcloudfrac(ngridmx)
    114 
    115 
    116 c=======================================================================
    117114
    118115c=======================================================================
    … …  
    613610
    614611      DO idt=1,ndt
    615         IF (idt.eq.ndt-day_step-1) then       !test
     612        IF (idt.eq.ndt) then       !test
    616613         lastcall=.true.
    617614         call stellarlong(day*1.0,zls)

  • trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F

    r253 r486  
    3131      real*8 FZEROI(L_NSPECTI)
    3232      real*8 taugsurf(L_NSPECTI,L_NGAUSS-1), fzero
    33 
    34 !      real*8 BSURFtest ! by RW for test
    3533
    3634      real*8 fup_tmp(L_NSPECTI),fdn_tmp(L_NSPECTI)
    … …  
    7371!      NTS   = TSURF*10.0D0-499
    7472!      NTT   = TTOP *10.0D0-499
    75 
    76 !      BSURFtest=0.0
    7773
    7874      DO 501 NW=1,L_NSPECTI
Note: See TracChangeset for help on using the changeset viewer.