Changeset 2599 for trunk/LMDZ.MARS


Ignore:
Timestamp:
Jan 4, 2022, 9:04:10 AM (3 years ago)
Author:
cmathe
Message:

MARS: clean co2condens_mod.F and remove dqsurf duplication after call co2condens

Location:
trunk/LMDZ.MARS
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r2597 r2599  
    35563556Also cleaned up and commented comsaison_h in the process.
    35573557
     3558== 04/01/2021 == CM
     3559- Clean co2condens_mod.F
     3560- remove dqsurf duplication after call co2condens

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

    r2566 r2599  
    323323      ENDDO
    324324     
     325      condens_layer(:,:) = zcondicea(:,:)
     326      condens_column(:) = 0.
     327
    325328      if (scavco2cond) then
    326329        call scavenging_by_co2(ngrid,nlayer,nq,ptimestep,pplev,zq,
    … …  
    339342        DO ig=1,ngrid
    340343          condens_column(ig) = sum(condens_layer(ig,:))
     344          zfallice(ig,1) = zdqssed_co2(ig)
    341345          piceco2(ig) = piceco2(ig) + zdqssed_co2(ig)*ptimestep
    342346        ENDDO
    … …  
    347351     &         " ",3,pdtc)
    348352     
    349        call WRITEdiagfi(ngrid,"zcondicea",
     353       call WRITEdiagfi(ngrid,"condens_layer",
    350354     &         "",
    351      &         " ",3,zcondicea)         
     355     &         " ",3,condens_layer)
    352356     
    353357       call WRITEdiagfi(ngrid,"zfallice",
    … …  
    374378c
    375379c     No microphysics of CO2 clouds
    376       if (.not.co2clouds)then
    377380      DO ig=1,ngrid
    378381         IF(latitude(ig).lt.0) THEN
    … …  
    383386            icap=1
    384387         ENDIF
    385        
    386388c
    387389c        Loop on where we have  condensation/ sublimation
    … …  
    390392     $      ((ztsrf(ig) .GT. ztcondsol(ig)) .AND.  ! ground sublim.
    391393     $      ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
    392             condsub(ig) = .true.
    393 
    394             IF (zfallice(ig,1).gt.0) then 
    395                  zfallheat=zfallice(ig,1)*
    396      &           (pphi(ig,1)- phisfi(ig) +
    397      &           cpice*(ztcond(ig,1)-ztcondsol(ig)))/latcond
     394            condsub(ig) = .true.
     395
     396            IF (zfallice(ig,1).gt.0 .and. .not. co2clouds) then
     397              zfallheat = zfallice(ig,1) * (pphi(ig,1)- phisfi(ig) +
     398     &                    cpice*(ztcond(ig,1)-ztcondsol(ig)))/latcond
    398399            ELSE
    399                  zfallheat=0.
     400              zfallheat = 0.
    400401            ENDIF
    401402c           condensation or partial sublimation of CO2 ice
    402403c           """""""""""""""""""""""""""""""""""""""""""""""
    403             zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig))
    404      &      /(latcond*ptimestep)         - zfallheat
     404            zcondices(ig) = pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig))
     405     &                      /(latcond*ptimestep)  - zfallheat
    405406            pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
    406             zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
    407407#ifdef MESOSCALE
    408408      print*, "not enough CO2 tracer in 1st layer to condense"
    … …  
    413413c           """"""""""""""""""""""""""""""""""""""""""""""""""""""
    414414            IF(ico2.ne.0) then
    415 c              Available CO2 tracer in layer 1 at end of timestep (kg/m2)
     415c             Available CO2 tracer in layer 1 at end of timestep (kg/m2)
    416416#ifndef MESOSCALE
    417                 availco2= pq(ig,1,ico2)*((ap(1)-ap(2))+
    418      &          (bp(1)-bp(2))*(pplev(ig,1)/g-zdiceco2(ig)*ptimestep))
     417                availco2 = pq(ig,1,ico2)*((ap(1)-ap(2))+
     418     &                     (bp(1)-bp(2))*(pplev(ig,1)/g -
     419     &                     (zcondices(ig) + zfallice(ig,1))*ptimestep))
     420                if ((zcondices(ig) + condens_layer(ig,1))*ptimestep
     421     &           .gt.availco2) then
     422                  zcondices(ig) = availco2/ptimestep -
     423     &                            condens_layer(ig,1)
     424                  pdtsrfc(ig) = (latcond/pcapcal(ig))*
     425     &                          (zcondices(ig)+zfallheat)
     426                end if
    419427#else
    420428                availco2 = pq(ig,1,igcm_co2)
    … …  
    422430     &                  ' CORRECTED FROM SIGMA LEVELS"
    423431#endif
    424                IF ((zcondices(ig) + zcondicea(ig,1))*ptimestep
    425      &           .gt.availco2) then
    426                    zcondices(ig) = availco2/ptimestep -zcondicea(ig,1)
    427                    zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
    428                    pdtsrfc(ig)=(latcond/pcapcal(ig))*
    429      &                          (zcondices(ig)+zfallheat)
    430                ENDIF
    431432            ENDIF
    432433#endif
    … …  
    437438            IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
    438439     &          -zcondices(ig))THEN
    439                zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1) 
     440               zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1)
    440441               pdtsrfc(ig)=(latcond/pcapcal(ig))*
    441442     &         (zcondices(ig)+zfallheat)
    442                zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
    443443            END IF
    444444
    445445c           Changing CO2 ice amount and pressure :
    446446c           """"""""""""""""""""""""""""""""""""
     447            zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
     448     &        + condens_column(ig)
     449            if (co2clouds) then
     450             ! add here only direct condensation/sublimation
     451            piceco2(ig) = piceco2(ig) + zcondices(ig)*ptimestep
     452            else
     453             ! add here also CO2 ice in the atmosphere
    447454            piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
     455            end if
    448456            pdpsrf(ig) = -zdiceco2(ig)*g
    449457
    … …  
    461469         END IF ! if there is condensation/sublimation
    462470      ENDDO ! of DO ig=1,ngrid
    463       else
    464         do ig=1,ngrid
    465           if ((ztsrf(ig)<ztcondsol(ig)) .OR. (zdqssed_co2(ig)/=0.) .OR.
    466      &       ((ztsrf(ig)>ztcondsol(ig)) .AND. (piceco2(ig)/=0.))) then
    467              condsub(ig) = .true.
    468              zfallheat = 0.
    469              !----------------------------------------------------------------------------------------------------------------------!
    470              ! Compute direct condensation/sublimation of CO2 ice
    471              !----------------------------------------------------------------------------------------------------------------------!
    472              zcondices(ig) = pcapcal(ig) * (ztcondsol(ig)-ztsrf(ig)) /
    473      &                       (latcond*ptimestep) - zfallheat
    474              pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig)) / ptimestep
    475              !----------------------------------------------------------------------------------------------------------------------!
    476              ! If there is not enough CO2 tracer in 1st layer to condense
    477              !----------------------------------------------------------------------------------------------------------------------!
    478              !     Available CO2 tracer in layer 1 at end of timestep (kg/m2)
    479 #ifndef MESOSCALE
    480              availco2 = pq(ig,1,igcm_co2) * ( (ap(1)-ap(2)) +
    481      &                  (bp(1)-bp(2)) * (pplev(ig,1)/g - zcondices(ig)
    482      &                  *ptimestep) )
    483 #else
    484              availco2 = pq(ig,1,igcm_co2)
    485              PRINT*, "MESOSCALE: CO2 tracer AT FIRST LEVEL IS NOT CORRECTED FROM SIGMA LEVELS"
    486 #endif
    487              if ( zcondices(ig) * ptimestep>availco2 ) then
    488                zcondices(ig) = availco2/ptimestep
    489                zdiceco2(ig) = zcondices(ig) + zdqssed_co2(ig)
    490                pdtsrfc(ig) = (latcond/pcapcal(ig)) *
    491      &                       (zcondices(ig)+zfallheat)
    492              end if
    493              !----------------------------------------------------------------------------------------------------------------------!
    494              ! If the entire CO2 ice layer sublimes (including what has just condensed in the atmosphere)
    495              !----------------------------------------------------------------------------------------------------------------------!
    496              if ( (piceco2(ig)/ptimestep) <= -zcondices(ig) ) then
    497                 zcondices(ig) = -piceco2(ig)/ptimestep
    498                 pdtsrfc(ig) = (latcond/pcapcal(ig)) *
    499      &                        (zcondices(ig)+zfallheat)
    500              end if
    501              !----------------------------------------------------------------------------------------------------------------------!
    502              ! Changing CO2 ice amount and pressure
    503              !----------------------------------------------------------------------------------------------------------------------!
    504              zdiceco2(ig) = zcondices(ig) + zdqssed_co2(ig) +
    505      &                      condens_column(ig)
    506              piceco2(ig) = piceco2(ig) + zcondices(ig)*ptimestep
    507              pdpsrf(ig) = -zdiceco2(ig) * g
    508 
    509              if (abs(pdpsrf(ig)*ptimestep)>pplev(ig,1)) then
    510              print *, 'STOP in condens'
    511              print *, 'condensing more than total mass'
    512              print *, 'Grid point ', ig
    513              print *, 'Longitude(degrees): ', longitude_deg(ig)
    514              print *, 'Latitude(degrees): ', latitude_deg(ig)
    515              print *, 'Ps = ', pplev(ig,1)
    516              print *, 'd Ps = ', pdpsrf(ig)
    517              call abort_physic('co2condens',
    518      &                         'condensing more than total mass', 1)
    519              end if
    520           end if
    521         end do
    522       end if
     471
    523472c  ********************************************************************
    524473c  Surface albedo and emissivity of the surface below the snow (emisref)
    … …  
    565514c  Mass fluxes through the sigma levels (kg.m-2.s-1)  (>0 when up)
    566515c  """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
    567             if (.not. co2clouds) then
    568               zmflux(1) = -zcondices(ig)
     516              zmflux(1) = -zcondices(ig) - zdqssed_co2(ig)
    569517              DO l=1,nlayer
    570                 zmflux(l+1) = zmflux(l) -zcondicea(ig,l)
     518                zmflux(l+1) = zmflux(l) - condens_layer(ig,l)
    571519#ifndef MESOSCALE
    572      &          + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
     520     &          + (bp(l)-bp(l+1))*(-pdpsrf(ig)/g)
    573521c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
    574522                if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) then
    … …  
    576524                end if
    577525#else
    578                 zmflux(l+1) = zmflux(l) - zcondicea(ig,l)
     526                zmflux(l+1) = zmflux(l) - condens_layer(ig,l)
    579527                if (abs(zmflux(l+1)).lt.1E-13) zmflux(l+1)=0.
    580528                PRINT*, "MESOSCALE: FLUX THROUGH SIGMA LEVELS from"//
    … …  
    582530#endif
    583531              END DO
    584             else
    585               zmflux(1) = - zcondices(ig) - zdqssed_co2(ig)
    586               do l = 1, nlayer
    587 #ifndef MESOSCALE
    588                 zmflux(l+1) = zmflux(l) - condens_layer(ig,l)
    589      &                        + (bp(l)-bp(l+1)) * (-pdpsrf(ig)/g)
    590 ! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
    591                 if (abs(zmflux(l+1))<1E-13.OR.bp(l+1)==0.) then
    592                   zmflux(l+1) = 0.
    593                 end if
    594 #else
    595                 zmflux(l+1) = zmflux(l) - condens_layer(ig,l)
    596                 if (abs(zmflux(l+1))<1E-13) then
    597                   zmflux(l+1) = 0.
    598                 end if
    599                 PRINT*, "MESOSCALE: FLUX THROUGH SIGMA LEVELS from"//
    600      &            "dPS HAVE TO BE IMPLEMENTED"
    601 #endif
    602               end do
    603             end if
    604532#ifdef MESOSCALE
    605533         print*,"absurd mass set because bp not available"
    … …  
    705633     &        ( zmflux(l)*(ztm(l) - ztc(l))
    706634     &        - zmflux(l+1)*(ztm(l+1) - ztc(l))
    707      &        + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l))  )
     635     &        + condens_layer(ig,l)*(ztcond(ig,l)-ztc(l))  )
    708636               ELSE
    709637                zdtsig(ig,l) = (1/masse(l)) *
    … …  
    733661#endif
    734662
    735             if (.not. co2clouds) then
    736663              do iq=1,nq
    737664!                if (noms(iq).eq.'co2') then
    738                 if (iq.eq.ico2) then
    739 c                 SPECIAL Case when the tracer IS CO2 :
    740                   DO l=1,nlayer
    741                     pdqc(ig,l,iq)= (1/masse(l)) *
    742      &              ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
    743      &              - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
    744      &              + zcondicea(ig,l)*(zqc(l,iq)-1.) )
    745                   END DO
    746                 else
    747                   DO l=1,nlayer
    748                     pdqc(ig,l,iq)= (1/masse(l)) *
    749      &             ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
    750      &             - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
    751      &             + zcondicea(ig,l)*zqc(l,iq) )
    752 
    753                     pdqc(ig,l,iq)=pdqc(ig,l,iq)+zdq_scav(ig,l,iq) ! zdq_scav=0 if watercloud=false
    754                   END DO
    755                 end if
    756               enddo
    757             else
    758 c           Tendencies on Q
    759               do iq=1,nq
    760 !               if (noms(iq).eq.'co2') then
    761665                if (iq.eq.ico2) then
    762666c                 SPECIAL Case when the tracer IS CO2 :
    … …  
    770674                  DO l=1,nlayer
    771675                    pdqc(ig,l,iq)= (1/masse(l)) *
    772      &              ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
    773      &              - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
    774      &              + condens_layer(ig,l)*zqc(l,iq) )
     676     &             ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
     677     &             - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
     678     &             + condens_layer(ig,l)*zqc(l,iq) )
     679
     680                    pdqc(ig,l,iq)=pdqc(ig,l,iq)+zdq_scav(ig,l,iq) ! zdq_scav=0 if watercloud=false
    775681                  END DO
    776682                end if
    777683              enddo
    778             end if
    779 
    780684
    781685       end if ! if (condsub)
    … …  
    785689c   CO2 snow / clouds scheme
    786690c ***************************************************************
    787       if (.not. co2clouds) then
    788691        call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
    789      &               zcondicea,zcondices,zfallice,pemisurf)
    790       else
    791         call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
    792      &               zcondicea,zcondices,zdqssed_co2,pemisurf)
    793       end if
     692     &               condens_layer,zcondices,zfallice,pemisurf)
    794693c ***************************************************************
    795694c Ecriture des diagnostiques

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

    r2597 r2599  
    874874          co2totA = co2totA + co2ice(ig)
    875875        end do
     876      else
     877        co2totA = 0.
     878        do ig=1,ngrid
     879          do l=1,nlayer
     880             co2totA = co2totA + (zplev(ig,l)-zplev(ig,l+1))/g*
     881     &             (pq(ig,l,igcm_co2)
     882     &        +pdq(ig,l,igcm_co2)*ptimestep)
     883          end do
     884          co2totA = co2totA + co2ice(ig)
     885        end do
    876886      endif ! of if (igcm_co2_ice.ne.0)
    877887c-----------------------------------------------------------------------
    … …  
    20592069            ENDDO
    20602070           ENDDO
    2061 
    2062           DO iq=1, nq
    2063             DO ig=1,ngrid
    2064               dqsurf(ig,iq)=dqsurf(ig,iq)+zdqsc(ig,iq)
    2065             ENDDO  ! (ig)
    2066            ENDDO    ! (iq)
    2067 
    20682071         ENDIF ! of IF (tracer)
    20692072
    … …  
    39533956          co2totB = co2totB + co2ice(ig)
    39543957        enddo
    3955 
     3958      else
     3959        co2totB = 0. ! added by C.M.
     3960        do ig=1,ngrid
     3961          do l=1,nlayer
     3962            co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g*
     3963     &             (pq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2)*ptimestep)
     3964          enddo
     3965          co2totB = co2totB + co2ice(ig)
     3966        enddo
     3967      endif ! of if (igcm_co2_ice.ne.0)
    39563968        call WRITEDIAGFI(ngrid,'co2conservation',
    39573969     &                     'Total CO2 mass conservation in physic',
    39583970     &                     '%',0,[(co2totA-co2totB)/co2totA])
    3959       endif ! of if (igcm_co2_ice.ne.0)
    39603971! XIOS outputs
    39613972#ifdef CPP_XIOS     
Note: See TracChangeset for help on using the changeset viewer.