Ignore:
Timestamp:
Mar 24, 2026, 5:27:57 PM (7 weeks ago)
Author:
yluo
Message:

Mars PCM:
Refactored the tracer mass-conservation fixer for dynamics introduced in r4142:
Added a runtime flag to enable/disable the correction of tracer mass non-conservation in dynamics.
Reorganized the implementation into a modular structure.
Restricted the correction to a subset of tracers: CO, O2, H2, HO2, H2O2, N2, Ar, and He.
Enabled correction for dynamics timesteps preceding the first physics timestep of a simulation.
YCL

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
1 added
6 edited

Legend:

Unmodified
Added
Removed

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

    r4102 r4150  
    5050     &      coal_coeff_mode,full_coag_equations
    5151      use lwxn_mod, only: lwxn_linear, lwxn_alphan, lwxn_ncouche
     52      use tracer_mass_fixer_dyn_mod, only: call_mass_fixer_dyn
    5253      use callkeys_mod, only: startphy_file, activice, activeco2ice,
    5354     &                        calladj, callatke, callcond,
    … …  
    12581259        endif
    12591260
     1261c        ----------------------------------------------------------
     1262
     1263! Tracer mass fixer for dynamics
     1264
     1265         write(*,*) "call tracer mass fixer for dynamics?"
     1266         call getin_p("call_mass_fixer_dyn", call_mass_fixer_dyn)
     1267         write(*,*) "call_mass_fixer_dyn = ", call_mass_fixer_dyn
     1268
    12601269! Test of incompatibility:
    12611270! if photochem is used, then water should be used too

  • trunk/LMDZ.MARS/libf/phymars/iostart.F90

    r3702 r4150  
    88   
    99    ! restartfi.nc file dimension identifiers: (see open_restartphy())
     10    INTEGER,SAVE :: idim0  ! "scalar" dimension
    1011    INTEGER,SAVE :: idim1  ! "index" dimension
    1112    INTEGER,SAVE :: idim2  ! "physical_points" dimension
    … …  
    510511        write(*,*)'open_restartphy: problem writing title '
    511512        write(*,*)trim(nf90_strerror(ierr))
     513      ENDIF
     514
     515      ierr=NF90_DEF_DIM(nid_restart,"scalar",1,idim0)
     516      IF (ierr/=NF90_NOERR) THEN
     517        write(*,*)'open_restartphy: problem defining scalar dimension '
     518        write(*,*)trim(nf90_strerror(ierr))
     519        CALL abort_physic("open_restartphy","Failed defining scalar",1)
    512520      ENDIF
    513521
    … …  
    10611069        ENDIF
    10621070        return ! nothing left to do
     1071      ELSEIF (var_size==1) THEN
     1072        ! We know it is a scalar
     1073        idim1d=idim0
    10631074      ELSEIF (var_size==length) THEN
    10641075        ! We know it is a "controle" kind of 1D array

  • trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90

    r3974 r4150  
    2626use nonoro_gwd_mix_mod,  only: du_eddymix_gwd, dv_eddymix_gwd, de_eddymix_rto, &
    2727                               df_eddymix_flx, dh_eddymix_gwd, dq_eddymix_gwd
     28use tracer_mass_fixer_dyn_mod, only: mass_predyn_co, mass_predyn_o2, mass_predyn_h2, mass_predyn_ho2,  &
     29                                     mass_predyn_h2o2, mass_predyn_n2, mass_predyn_ar, mass_predyn_he, &
     30                                     found_startfi_co, found_startfi_o2, found_startfi_h2, found_startfi_ho2, &
     31                                     found_startfi_h2o2, found_startfi_n2, found_startfi_ar, found_startfi_he, &
     32                                     call_mass_fixer_dyn
    2833use compute_dtau_mod,    only: dtau
    2934use dust_rad_adjust_mod, only: dust_rad_adjust_prev,dust_rad_adjust_next
    … …  
    718723             minval(df_eddymix_flx ), maxval(df_eddymix_flx )
    719724
     725! Correction to tracer mass non-conservation in dynamics
     726if (startphy_file .and. call_mass_fixer_dyn) then
     727   call get_var("mass_predyn_co", mass_predyn_co, found)
     728   if (found) then
     729      found_startfi_co = .true.
     730   else
     731      write(*,*) "phyetat0: <mass_predyn_co> not in file"
     732   endif
     733
     734   call get_var("mass_predyn_o2", mass_predyn_o2, found)
     735   if (found) then
     736      found_startfi_o2 = .true.
     737   else
     738      write(*,*) "phyetat0: <mass_predyn_o2> not in file"
     739   endif
     740
     741   call get_var("mass_predyn_h2", mass_predyn_h2, found)
     742   if (found) then
     743      found_startfi_h2 = .true.
     744   else
     745      write(*,*) "phyetat0: <mass_predyn_h2> not in file"
     746   endif
     747
     748   call get_var("mass_predyn_ho2", mass_predyn_ho2, found)
     749   if (found) then
     750      found_startfi_ho2 = .true.
     751   else
     752      write(*,*) "phyetat0: <mass_predyn_ho2> not in file"
     753   endif
     754
     755   call get_var("mass_predyn_h2o2", mass_predyn_h2o2, found)
     756   if (found) then
     757      found_startfi_h2o2 = .true.
     758   else
     759      write(*,*) "phyetat0: <mass_predyn_h2o2> not in file"
     760   endif
     761
     762   call get_var("mass_predyn_n2", mass_predyn_n2, found)
     763   if (found) then
     764      found_startfi_n2 = .true.
     765   else
     766      write(*,*) "phyetat0: <mass_predyn_n2> not in file"
     767   endif
     768
     769   call get_var("mass_predyn_ar", mass_predyn_ar, found)
     770   if (found) then
     771      found_startfi_ar = .true.
     772   else
     773      write(*,*) "phyetat0: <mass_predyn_ar> not in file"
     774   endif
     775
     776   call get_var("mass_predyn_he", mass_predyn_he, found)
     777   if (found) then
     778      found_startfi_he = .true.
     779   else
     780      write(*,*) "phyetat0: <mass_predyn_he> not in file"
     781   endif
     782endif ! if (startphy_file .and. call_mass_fixer_dyn)
     783
    720784! tracer on surface
    721785if (nq.ge.1) then

  • trunk/LMDZ.MARS/libf/phymars/phyredem.F90

    r3964 r4150  
    198198  use iostart,             only: open_restartphy, close_restartphy, put_var, put_field
    199199  use tracer_mod,          only: noms ! tracer names
     200  use tracer_mod,          only: igcm_co, igcm_o2, igcm_h2, igcm_ho2, igcm_h2o2, igcm_n2, igcm_ar, igcm_he
    200201  use nonoro_gwd_ran_mod,  only: du_nonoro_gwd, dv_nonoro_gwd
    201202  use compute_dtau_mod,    only: dtau
    … …  
    206207  use comslope_mod,        only: nslope
    207208  use paleoclimate_mod,    only: paleoclimate
     209  use tracer_mass_fixer_dyn_mod, only: call_mass_fixer_dyn, mass_predyn
    208210  use callkeys_mod,        only: calltherm, dustinjection, calllott_nonoro
    209211  use callkeys_mod, only: CLFvarying
    … …  
    392394  endif
    393395
     396  ! Correction to tracer mass non-conservations in dynamics
     397  if (call_mass_fixer_dyn) then
     398    if (igcm_co .ne. 0) then
     399      call put_var("mass_predyn_co", "Global mass of CO before the last physiq step ends", mass_predyn(igcm_co))
     400    endif
     401    if (igcm_o2 .ne. 0) then
     402      call put_var("mass_predyn_o2", "Global mass of O2 before the last physiq step ends", mass_predyn(igcm_o2))
     403    endif
     404    if (igcm_h2 .ne. 0) then
     405      call put_var("mass_predyn_h2", "Global mass of H2 before the last physiq step ends", mass_predyn(igcm_h2))
     406    endif
     407    if (igcm_ho2 .ne. 0) then
     408      call put_var("mass_predyn_ho2", "Global mass of HO2 before the last physiq step ends", mass_predyn(igcm_ho2))
     409    endif
     410    if (igcm_h2o2 .ne. 0) then
     411      call put_var("mass_predyn_h2o2", "Global mass of H2O2 before the last physiq step ends", mass_predyn(igcm_h2o2))
     412    endif
     413    if (igcm_n2 .ne. 0) then
     414      call put_var("mass_predyn_n2", "Global mass of N2 before the last physiq step ends", mass_predyn(igcm_n2))
     415    endif
     416    if (igcm_ar .ne. 0) then
     417      call put_var("mass_predyn_ar", "Global mass of Ar before the last physiq step ends", mass_predyn(igcm_ar))
     418    endif
     419    if (igcm_he .ne. 0) then
     420      call put_var("mass_predyn_he", "Global mass of He before the last physiq step ends", mass_predyn(igcm_he))
     421    endif
     422  endif
     423
    394424  ! Close file
    395425  call close_restartphy

  • trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90

    r4009 r4150  
    6565      use comslope_mod, ONLY: nslope,end_comslope_h,ini_comslope_h
    6666      use paleoclimate_mod, ONLY: end_paleoclimate_h,ini_paleoclimate_h
     67      use tracer_mass_fixer_dyn_mod, only: ini_tracer_mass_fixer_dyn, &
     68                                           end_tracer_mass_fixer_dyn
    6769      use netcdf
    6870#ifndef MESOSCALE
    … …  
    200202      call ini_paleoclimate_h(ngrid,nslope)
    201203
     204      ! allocate arrays in "tracer_mass_fixer_dyn_mod"
     205      call end_tracer_mass_fixer_dyn
     206      call ini_tracer_mass_fixer_dyn(nq)
     207
    202208      END SUBROUTINE phys_state_var_init
    203209

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

    r4143 r4150  
    275275      REAL,INTENT(in) :: pt(ngrid,nlayer) ! temperature (K)
    276276      REAL,INTENT(in) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
    277       REAL pq_corrdyn(ngrid,nlayer,nq) ! tracer mass mixing ratio after correcting
    278                                        ! mass non-conservations in dynamics
    279277      REAL,INTENT(in) :: flxw(ngrid,nlayer) ! vertical mass flux (ks/s) at lower boundary of layer
    280278
    … …  
    286284      REAL,INTENT(out) :: pdt(ngrid,nlayer) ! temperature tendency (K/s)
    287285      REAL,INTENT(out) :: pdq(ngrid,nlayer,nq) ! tracer tendencies (../kg/s)
    288       REAL pdq_corrdyn(ngrid,nlayer,nq) ! tracer mass mixing ratio tendencies in order to correct
    289                                         ! mass non-conservations in dynamics
    290286      REAL,INTENT(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s)
    291287
    … …  
    451447      real colden(ngrid,nq)     ! vertical column of tracers
    452448      real mass(nq)             ! global mass of tracers (g)
    453       real, allocatable, save :: mass_predyn(:) ! tracer mass in the whole atmosphere
    454                                                 ! before dynamics
    455       real mass_postdyn(nq) ! tracer mass in the whole atmosphere after dynamics
    456 
    457 !$OMP THREADPRIVATE(mass_predyn)
    458 
    459       real masscorrfac(nq) ! mass correction factor for non-conservations in dynamics
    460449      REAL mtot(ngrid)          ! Total mass of water vapor (kg/m2)
    461450      REAL mstormdtot(ngrid)    ! Total mass of stormdust tracer (kg/m2)
    … …  
    826815     &                              presnivs,pseudoalt,mlayer)
    827816#endif
    828 
    829          allocate(mass_predyn(nq))
    830 
    831817      ENDIF        !  (end of "if firstcall")
    832818
    … …  
    892878      zplay(:,:) = pplay(:,:)
    893879      ps(:) = pplev(:,1)
    894 
    895 c     Correct mass non-conservation in dynamics
    896 c     ~~~~~~~~~~~~~~~~~~
    897       call compute_global_mass(pq, ngrid, nlayer, nq, zplev,
    898      &                         mass_postdyn)
    899 
    900       if (.not.firstcall) then ! At first call, mass_predyn does not exist
    901          do iq = 1, nq
    902             if (noms(iq) .ne. "dust_mass" .and.
    903      &          noms(iq) .ne. "dust_number" .and.
    904      &          noms(iq) .ne. "ccn_mass" .and.
    905      &          noms(iq) .ne. "ccn_number" .and.
    906      &          noms(iq) .ne. "ccnco2_mass" .and.
    907      &          noms(iq) .ne. "ccnco2_number" .and.
    908      &          noms(iq) .ne. "stormdust_mass" .and.
    909      &          noms(iq) .ne. "stormdust_number" .and.
    910      &          noms(iq) .ne. "topdust_mass" .and.
    911      &          noms(iq) .ne. "topdust_number") then
    912                masscorrfac(iq) = mass_predyn(iq)/mass_postdyn(iq)
    913             endif
    914          enddo
    915 
    916          pq_corrdyn = pq
    917          do ig = 1, ngrid
    918             do l = 1, min(nlayer, 40) ! Correction only applied to the lowest 40 layers
    919                do iq = 1, nq
    920                   if (noms(iq) .ne. "dust_mass" .and.
    921      &                noms(iq) .ne. "dust_number" .and.
    922      &                noms(iq) .ne. "ccn_mass" .and.
    923      &                noms(iq) .ne. "ccn_number" .and.
    924      &                noms(iq) .ne. "ccnco2_mass" .and.
    925      &                noms(iq) .ne. "ccnco2_number" .and.
    926      &                noms(iq) .ne. "stormdust_mass" .and.
    927      &                noms(iq) .ne. "stormdust_number" .and.
    928      &                noms(iq) .ne. "topdust_mass" .and.
    929      &                noms(iq) .ne. "topdust_number" .and.
    930      &                iq .ne. igcm_co2) then
    931                      ! Correcting species other than CO2
    932                      pq_corrdyn(ig, l, iq) =
    933      &                    pq(ig, l, iq) * masscorrfac(iq)
    934                      ! Adjust CO2 to compensate for mass creation/loss
    935                      pq_corrdyn(ig, l, igcm_co2) =
    936      &                    pq_corrdyn(ig, l, igcm_co2) -
    937      &                    pq(ig, l, iq) * (masscorrfac(iq) - 1.)
    938                   endif
    939                enddo
    940             enddo
    941          enddo
    942 
    943          pdq_corrdyn = (pq_corrdyn - pq)/ptimestep ! Mass fixer pdq
    944          pdq = pdq + pdq_corrdyn
    945       endif
    946880
    947881!#ifndef MESOSCALE
    … …  
    26342568        ENDDO
    26352569      ENDDO
    2636 
    2637       call compute_global_mass(zq, ngrid, nlayer, nq, zplev,
    2638      &                         mass_predyn)
    26392570
    26402571      ! Density
    … …  
    41994130      END SUBROUTINE physiq
    42004131
    4201       subroutine compute_global_mass(zq, ngrid, nlayer, nq, zplev, mass)
    4202 
    4203          use tracer_mod, only: mmol
    4204          use comcstfi_h, only: g
    4205          use geometry_mod, only: cell_area
    4206          use tracer_mod, only: noms
    4207          use planetwide_mod, only: planetwide_sumval
    4208 
    4209          real, intent(in)    :: zq(ngrid, nlayer, nq)
    4210          integer, intent(in) :: ngrid
    4211          integer, intent(in) :: nlayer
    4212          integer, intent(in) :: nq
    4213          real, intent(in)    :: zplev(ngrid, nlayer + 1)
    4214          real, intent(out)   :: mass(nq)
    4215          real    colden(ngrid, nq)
    4216          integer iq, ig, l
    4217 
    4218          do iq = 1, nq
    4219             if (noms(iq) .ne. "dust_mass" .and.
    4220      &          noms(iq) .ne. "dust_number" .and.
    4221      &          noms(iq) .ne. "ccn_mass" .and.
    4222      &          noms(iq) .ne. "ccn_number" .and.
    4223      &          noms(iq) .ne. "ccnco2_mass" .and.
    4224      &          noms(iq) .ne. "ccnco2_number" .and.
    4225      &          noms(iq) .ne. "stormdust_mass" .and.
    4226      &          noms(iq) .ne. "stormdust_number" .and.
    4227      &          noms(iq) .ne. "topdust_mass" .and.
    4228      &          noms(iq) .ne. "topdust_number") then
    4229 
    4230                do ig = 1, ngrid
    4231                    colden(ig, iq) = 0.
    4232                end do
    4233                do l = 1, nlayer
    4234                   do ig = 1, ngrid
    4235                      colden(ig, iq) = colden(ig, iq) + zq(ig, l, iq)*
    4236      &                                (zplev(ig, l) - zplev(ig, l+1))*
    4237      &                                6.022e22/(mmol(iq)*g)
    4238                   end do
    4239                end do
    4240 
    4241                call planetwide_sumval(colden(:, iq)/6.022e23
    4242      &                    *mmol(iq)*1.e4*cell_area, mass(iq))
    4243 
    4244             end if
    4245          end do
    4246 
    4247       end subroutine
    4248 
    42494132      END MODULE physiq_mod
Note: See TracChangeset for help on using the changeset viewer.