Ignore:
Timestamp:
Oct 18, 2018, 12:39:43 PM (6 years ago)
Author:
emillour
Message:

Mars GCM:
Cleanup "newcondens" (remove obsolete options), change its name
to co2condens and turn it into a module.
EM

File:
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F

    r2007 r2009  
    1       SUBROUTINE newcondens(ngrid,nlayer,nq,ptimestep,
     1      MODULE co2condens_mod
     2
     3      IMPLICIT NONE
     4
     5      CONTAINS
     6
     7      SUBROUTINE co2condens(ngrid,nlayer,nq,ptimestep,
    28     $                  pcapcal,pplay,pplev,ptsrf,pt,
    39     $                  pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq,
    … …  
    1117       use surfdat_h, only: emissiv, phisfi
    1218       use geometry_mod, only: latitude ! grid point latitudes (rad)
    13        use planete_h
    14        USE comcstfi_h
     19       use planete_h, only: obliquit
     20       use comcstfi_h, only: cpp, g, r, pi
    1521#ifndef MESOSCALE
    1622       USE vertical_layers_mod, ONLY: bp
    … …  
    2430c   (Scheme described in Forget et al., Icarus, 1998)
    2531c
    26 c   author:   Francois Forget     1994-1996
     32c   author:   Francois Forget     1994-1996 ; updated 1996 -- 2018
    2733c   ------
    2834c             adapted to external CO2 ice clouds scheme by Deborah Bardet (2018) '
    29 c
    30 c   input:
    31 c   -----
    32 c    ngrid                 nombre de points de verticales
    33 c                          (toutes les boucles de la physique sont au
    34 c                          moins vectorisees sur ngrid)
    35 c    nlayer                nombre de couches
    36 c    pplay(ngrid,nlayer)   Pressure levels
    37 c    pplev(ngrid,nlayer+1) Pressure levels
    38 c    pt(ngrid,nlayer)      temperature (en K)
    39 c    ptsrf(ngrid)          temperature de surface
    40 c
    41 c                    \
    42 c    pdt(ngrid,nlayer)\  derivee temporelle physique avant condensation
    43 c                     /  ou sublimation pour  pt,ptsrf
    44 c    pdtsrf(ngrid)   /
    45 c
    46 c   output:
    47 c   -------
    48 c
    49 c    pdpsrf(ngrid)      \  derivee temporelle physique (contribution de
    50 c    pdtc(ngrid,nlayer) /  la condensation ou sublimation) pour Ps,pt,ptsrf
    51 c    pdtsrfc(ngrid)    /
    52 c
    53 c   Entree/sortie
    54 c   -------------
    55 c   
    56 c    piceco2(ngrid) :      quantite de glace co2 au sol (kg/m2)
    57 c    psolaralb(ngrid,2) :  albedo au sol
    58 c    pemisurf(ngrid)     :  emissivite du sol             
    59 
    6035c
    6136c=======================================================================
    … …  
    11893c    -----------------
    11994
    120 c   variables used for albedo parametrization       
    121 c --------------------------------------------     
    12295      INTEGER i,j
    123 c      REAL Fluxmean(jjp1)
    124       INTEGER l,ig,iq,icap,nmix
    125       LOGICAL transparency, fluxdependent
    126 c   flag transparency if you want to make the co2ice semi-transparent       
    127       PARAMETER(transparency=.true.)
    128 c   flag fluxdependent if you want the co2ice albedo to be dependent on
    129 c   the incident solar flux     
    130       PARAMETER(fluxdependent=.false.)
    131       REAL slopy,alpha,constA,constB,constT,albediceF_new(ngrid)
     96      INTEGER l,ig,iq,icap
    13297      REAL zt(ngrid,nlayer)
    13398      REAL zcpi
    … …  
    135100      REAL ztcondsol(ngrid) ! CO2 condensation temperature (surface)
    136101      REAL zdiceco2(ngrid)
    137       REAL zcondicea(ngrid,nlayer)
    138       REAL zcondices(ngrid)
    139       REAL zfallice(ngrid,nlayer+1) , zfallheat
     102      REAL zcondicea(ngrid,nlayer) ! condensation rate in layer  l  (kg/m2/s)
     103      REAL zcondices(ngrid) ! condensation rate on the ground (kg/m2/s)
     104      REAL zfallice(ngrid,nlayer+1) ! amount of ice falling from layer l (kg/m2/s)
     105      REAL zfallheat
    140106      REAL zmflux(nlayer+1)
    141107      REAL zu(nlayer),zv(nlayer)
    … …  
    161127! then condensation temperature is computed using partial pressure of CO2
    162128      logical,parameter :: improved_ztcond=.true.
    163 ! Bound co2 (tracer) values...
    164       logical,parameter :: bound_qco2=.false.
    165       real,parameter :: qco2max=1.1
    166       real,parameter :: qco2mini=0.1
    167       real :: zqco2
    168129
    169130c   local saved variables
    … …  
    174135      real,parameter :: cpice=1000. ! (J.kg-1.K-1) specific heat of CO2 ice
    175136      REAL,SAVE :: acond,bcond,ccond
    176 !      REAL,SAVE :: albediceF(ngrid)
    177137      real,save :: m_co2, m_noco2, A , B
    178138
    179139      LOGICAL,SAVE :: firstcall = .true. !,firstcall2=.true.
    180 
    181       integer flag
    182140
    183141c D.BARDET: to debug
    … …  
    199157         firstcall=.false.
    200158         write(*,*) 'Newcondens: improved_ztcond=',improved_ztcond
    201          write(*,*) 'Newcondens: bound_qco2=',bound_qco2
    202159         PRINT*,'In newcondens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
    203160         PRINT*,'acond,bcond,ccond',acond,bcond,ccond
    … …  
    454411c  Surface albedo and emissivity of the surface below the snow (emisref)
    455412c  ********************************************************************
    456 c     Prepare the case where albedo varies with insolation:
    457 c     ----------------------------------------------------
    458 !      if (fluxdependent) then
    459 !
    460 c       Calcul du flux moyen (zonal mean)
    461 !        do j=1,jjp1
    462 !           Fluxmean(j)=0     
    463 !           do i=1,iim
    464 !               ig=1+(j-2)*iim +i
    465 !               if(j.eq.1) ig=1
    466 !               if(j.eq.jjp1) ig=ngrid
    467 !               Fluxmean(j)=Fluxmean(j)+fluxsurf_sw(ig,1)
    468 !     $            +fluxsurf_sw(ig,2)
    469 !           enddo
    470 !           Fluxmean(j)=Fluxmean(j)/float(iim)
    471 !        enddo
    472 !
    473 c       const A and B used to calculate the albedo which depends on solar flux
    474 c       albedice=constA+constB*Flux
    475 c       constT = time step to calculate the solar flux when flux decreases         
    476 !         constA=0.26
    477 c        constA=0.33
    478 c        constA=0.186
    479 !         constB=0.00187
    480 !         constT=10
    481 !      endif ! of if (fluxdependent)
    482413
    483414!     Check that amont of CO2 ice is not problematic
    … …  
    493424!     ----------------------------------------
    494425      CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
    495 
    496 c     Calcul de l'albedo
    497 c     ------------------
    498 !         do ig =1,ngrid
    499 !           IF(ig.GT.ngrid/2+1) THEN
    500 !              icap=2
    501 !           ELSE
    502 !              icap=1
    503 !           ENDIF
    504 !           IF(firstcall2) THEN
    505 !           albediceF(ig)=albedice(icap)
    506 !           ENDIF
    507 c       if there is still co2ice        ccccccccccccccccccccccc
    508 !           if (piceco2(ig).gt.0) then
    509 !              emisref(ig) = emisice(icap)
    510 
    511 c     if flux dependent albedo is used
    512 c     --------------------------------
    513 !             if (fluxdependent) then
    514 !                j=INT((ig-2)/iim)+2
    515 !                if(ig.eq.1) j=1
    516 !                if(ig.eq.ngrid) j=jjp1
    517 c                albediceF_new(ig)=MIN(constA+constB*Fluxmean(j),
    518 c     $          constA+constB*250)
    519 !                  albediceF_new(ig)=constA+constB*Fluxmean(j)
    520 !                if (albediceF(ig).gt.albediceF_new(ig)) then
    521 !                   albediceF(ig)=albediceF(ig)+ ptimestep/(daysec*
    522 !     $             constT)*(albediceF_new(ig)-albediceF(ig))
    523 !                else
    524 !                   albediceF(ig)=albediceF_new(ig)
    525 !                endif
    526 c               if part of the ice is transparent
    527 c               slopy = pente de la droite:  alpha = y*co2ice/1620
    528 c               pour une valeur superieur a une epaisseur de glace donnee
    529 c               ici, epaisseur limite = 10cm
    530 !                if (transparency) then
    531 !                slopy=1/(1620*0.10)
    532 !                alpha=MIN(slopy*piceco2(ig),1.)
    533 !                psolaralb(ig,1) = alpha*albediceF(ig)
    534 !     $                           +(1-alpha)*albedodat(ig)
    535 !                psolaralb(ig,2) = psolaralb(ig,1)
    536 !                else
    537 !                psolaralb(ig,1) = albediceF(ig)
    538 !                psolaralb(ig,2) = psolaralb(ig,1)
    539 !                endif
    540 !             else
    541 c           transparency set to true and fluxdependent set to false           
    542 !                if (transparency) then
    543 !                slopy=1/(1620*0.10)
    544 !                alpha=MIN(slopy*piceco2(ig),1.)
    545 !                psolaralb(ig,1) = alpha*albedice(icap)
    546 !     $                           +(1-alpha)*albedodat(ig)
    547 !                psolaralb(ig,2) = psolaralb(ig,1)
    548 !                else
    549 c           simplest case: transparency and flux dependent set to false
    550 !                psolaralb(ig,1) = albedice(icap)
    551 !                psolaralb(ig,2) = albedice(icap)
    552 !                endif
    553 !             endif   
    554 c         no more co2ice, albedo = ground albedo
    555 !           else
    556 !              psolaralb(ig,1) = albedodat(ig)
    557 !              psolaralb(ig,2) = albedodat(ig)
    558 !              emisref(ig) = emissiv
    559 !              pemisurf(ig) = emissiv
    560 !           endif
    561 !         end do ! end of the ig loop
    562426
    563427! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
    … …  
    722586            enddo
    723587
    724 c           --------------------------------------------------------
    725 c           Roughly Simulate Molecular mixing when CO2 is too depleted by
    726 c           Surface condensation (mixing starts if qco2 < qco2min ) 
    727 c           FF 06/2004
    728 c           WARNING : this is now done in convadj, better (FF 02/2005)
    729 c           --------------------------------------------------------
    730             flag=0  ! now done in convadj : must be =0
    731             if (flag.eq.1) then
    732             if(ico2.gt.0) then   ! relevant only if one tracer is CO2
    733              if(pq(ig,1,ico2)+(pdq(ig,1,ico2)+pdqc(ig,1,ico2))*ptimestep
    734      &       .lt.qco2min) then
    735                 do iq=1,nq
    736                   zq(1,iq)=pq(ig,1,iq)
    737      &                     +(pdq(ig,1,iq)+pdqc(ig,1,iq))*ptimestep
    738                   Smq(1,iq) = masse(1)*zq(1,iq)
    739                 end do
    740                 Sm(1)  = masse(1)
    741                 do l =2,nlayer
    742                   do iq=1,nq
    743                      zq(l,iq)=pq(ig,l,iq)
    744      &                        +(pdq(ig,l,iq)+pdqc(ig,l,iq))*ptimestep
    745                      smq(l,iq) = smq(l-1,iq) + masse(l)*zq(l,iq)
    746                   end do
    747                   sm(l) = sm(l-1) + masse(l)
    748                   if(zq(l,ico2).gt.qco2min) then
    749 c                   mixmas: mass of atmosphere that must be mixed to reach qco2min
    750                     mixmas = (sm(l-1)*zq(l,ico2)-Smq(l-1,ico2))
    751      &                      /(zq(l,ico2)-qco2min)
    752                     if((mixmas.le.sm(l)))then
    753 c                      OK if mixed mass less than mass of layers affected
    754                        nmix=l   ! number of layer affected by mixing
    755                        goto 99
    756                     end if
    757                   end if
    758                  end do
    759  99              continue
    760                  do iq=1,nq
    761                    qmix=zq(nmix,iq)
    762      &             +(Smq(nmix-1,iq)-zq(nmix,iq)*Sm(nmix-1))/mixmas
    763                    do l=1,nmix-1
    764                       pdqc(ig,l,iq)=
    765      &                  (qmix-pq(ig,l,iq))/ptimestep - pdq(ig,l,iq)
    766                    end do
    767 c                  layer only partly mixed :
    768                    pdqc(ig,nmix,iq)=(
    769      &             qmix+(Sm(nmix)-mixmas)*(zq(nmix,iq)-qmix)/masse(nmix)
    770      &             -pq(ig,nmix,iq))/ptimestep -pdq(ig,nmix,iq)
    771 
    772                  end do
    773              end if
    774             end if
    775 
    776         endif   ! (flag.eq.1)
    777588       end if ! if (condsub)
    778589      END DO  ! loop on ig
    … …  
    808619!     &'Precipitation of co2 ice',
    809620!     & 'kg.m-2.s-1',2,zfallice(1,1))
    810 
    811 !! Specific stuff to bound co2 tracer values ....
    812       if (bound_qco2.and.(ico2.ne.0)) then
    813         do ig=1,ngrid
    814           do l=1,nlayer
    815             zqco2=pq(ig,l,ico2)
    816      &            +(pdq(ig,l,ico2)+pdqc(ig,l,ico2))*ptimestep
    817             if (zqco2.gt.qco2max) then
    818             ! correct pdqc:
    819               pdqc(ig,l,ico2)=((qco2max-pq(ig,l,ico2))/ptimestep)
    820      &                               -pdq(ig,l,ico2)
    821               write(*,*) "newcondens: adapting pdqc(ig,l,ico2)",
    822      &                   " so that co2 conc. does not exceed",qco2max
    823               write(*,*) "   ig:",ig,"  l:",l
    824             endif ! of if (zqco2.gt.qco2max)
    825             if (zqco2.lt.qco2mini) then
    826             ! correct pdqc:
    827               pdqc(ig,l,ico2)=((qco2mini-pq(ig,l,ico2))/ptimestep)
    828      &                               -pdq(ig,l,ico2)
    829               write(*,*) "newcondens: adapting pdqc(ig,l,ico2)",
    830      &                   " so that co2 conc. is not less than",qco2mini
    831               write(*,*) "   ig:",ig,"  l:",l
    832             endif ! of if (zqco2.lt.qco2mini)
    833           end do
    834         enddo
    835       endif ! of if (bound_qco2.and.(ico2.ne.0)) then
    836621
    837622#ifndef MESOSCALE
    … …  
    853638#endif
    854639
    855       end
     640      END SUBROUTINE co2condens
    856641
    857642
    … …  
    977762c         end if
    978763
    979       end
     764      END SUBROUTINE vl1d
     765     
     766      END MODULE co2condens_mod
Note: See TracChangeset for help on using the changeset viewer.