Changeset 4152 for trunk/LMDZ.COMMON


Ignore:
Timestamp:
Mar 25, 2026, 11:19:02 AM (6 weeks ago)
Author:
jbclement
Message:

PEM:

  • Fix outputs "diagevo.nc" for 3D data.
  • Fix sign in computing exchanges due to adsorption/ice table and in balancing H2O flux from/into atmosphere.
  • Few cleanings.

JBC

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/changelog.txt

    r4147 r4152  
    940940- Making the PEM able to do 0 year.
    941941- Explicit information about the frost values computed by the PEM + enforcing positivity of yearly minima.
     942
     943== 25/03/2026 == JBC
     944- Fix outputs "diagevo.nc" for 3D data.
     945- Fix sign in computing exchanges due to adsorption/ice table and in balancing H2O flux from/into atmosphere.
     946- Few cleanings.

  • trunk/LMDZ.COMMON/libf/evolution/clim_state_init.F90

    r4135 r4152  
    377377
    378378    ! H2O ice
    379     call print_msg("'h2o_ice' is initialized with default value 'h2oice_huge_ini' where 'watercaptag' is true and where yearly minimum of frost can be considered as a huge reservoir ('threshold_h2oice_cap').",LVL_NFO)
     379    call print_msg("'h2o_ice' is initialized with default value 'h2oice_huge_ini' where 'watercaptag' is true and where yearly minimum of frost can be considered as a huge reservoir ('threshold_h2oice_cap').",LVL_WRN)
    380380    h2o_ice(:,:) = 0._dp
    381381    do i = 1,ngrid
    … …  
    388388
    389389    ! CO2 ice
    390     call print_msg("'co2_ice' is initialized with 'perennial_co2ice' and yearly minimum of frost found in the PCM.",LVL_NFO)
     390    call print_msg("'co2_ice' is initialized with 'perennial_co2ice' and yearly minimum of frost found in the PCM.",LVL_WRN)
    391391    co2_ice(:,:) = co2_perice_PCM(:,:) + co2frost_PCM(:,:) - co2_frost4PCM(:,:)
    392392
    … …  
    486486    ! Layering
    487487    if (do_layering) then
    488         call print_msg('layerings_map is initialized with sub-surface strata.',LVL_NFO)
    489         call print_msg("Ice is added with 'h2oice_huge_ini' where 'watercaptag' is true and otherwise with 'perennial_co2ice' found in the PCM.",LVL_NFO)
     488        call print_msg('layerings_map is initialized with sub-surface strata.',LVL_WRN)
     489        call print_msg("Ice is added with 'h2oice_huge_ini' where 'watercaptag' is true and otherwise with 'perennial_co2ice' found in the PCM.",LVL_WRN)
    490490        do i = 1,ngrid
    491491            if (is_h2o_perice_PCM(i)) then

  • trunk/LMDZ.COMMON/libf/evolution/config.F90

    r4134 r4152  
    452452! Initialize physical data
    453453!     Arguments order: rad, g, mugaz, rcp
    454 call print_msg('    > Initializing physical constants',LVL_NFO)
     454call print_msg('> Initializing physical constants',LVL_NFO)
    455455call init_physics(controle(5),controle(7),controle(8),controle(9))
    456456
    457457! Initialize soil data
    458 call print_msg('    > Initializing soil parameters',LVL_NFO)
     458call print_msg('> Initializing soil parameters',LVL_NFO)
    459459volcapa = controle(35)
    460460if (abs(volcapa) < minieps) call stop_clean(__FILE__,__LINE__,'volcapa is 0 in "'//startfi_name//'"!',1)
    … …  
    462462! Initialize orbital data
    463463!     Arguments order: Obliquity, Perihelion, Aphelion, Date of perihelion, Year length, Sol length
    464 call print_msg('    > Initializing orbital characteristics of the planet',LVL_NFO)
     464call print_msg('> Initializing orbital characteristics of the planet',LVL_NFO)
    465465call init_orbit(controle(18),controle(15),controle(16),controle(17),controle(14),controle(10))
    466466

  • trunk/LMDZ.COMMON/libf/evolution/deftank/README

    r4110 r4152  
    6262The PEM simulation generates the following files:
    6363    > the usual outputs of the PCM: "restartfi.nc", "restart.nc"/"restart1D.txt", "diagfi.nc", etc;
    64     > the XIOS outputs of the PCM: "Xoutdaily4pem*.nc"/"Xoutyearly4pem*.nc";
     64    > the XIOS outputs of the PCM: "xoutdaily4pem*.nc"/"xoutyearly4pem*.nc";
    6565    > the outputs of the chained simulation: "pem_workflow.log", "pem_workflow.sts" and possibly "kill_pem_workflow.sh";
    6666    > the usual outputs of the PEM: "restartevo.nc", "restartfi.nc", "restart.nc"/"restart1D.txt" and "diagevo.nc".

  • trunk/LMDZ.COMMON/libf/evolution/deftank/clean.sh

    r4110 r4152  
    1919rm -f start*.nc
    2020rm -f used_*
    21 rm -f Xout*4pem*.nc
     21rm -f xout*4pem*.nc
    2222rm -f xios_client_*
    2323rm -f diag*.nc

  • trunk/LMDZ.COMMON/libf/evolution/deftank/pcm_run.job

    r4110 r4152  
    6565k=0
    6666if [ $(echo "$k > 0" | bc) -eq 1 ]; then # Only the last 2 years are taken for the PEM
    67     cp Xoutdaily4pem.nc Xoutdaily4pem_Y${k}.nc
    68     cp Xoutyearly4pem.nc Xoutyearly4pem_Y${k}.nc
     67    cp xoutdaily4pem.nc xoutdaily4pem_y${k}.nc
     68    cp xoutyearly4pem.nc xoutyearly4pem_y${k}.nc
    6969fi
    70 mv Xoutdaily4pem.nc diags/Xoutdaily4pem${i_pcm_run}.nc
    71 mv Xoutyearly4pem.nc diags/Xoutyearly4pem${i_pcm_run}.nc
     70mv xoutdaily4pem.nc diags/xoutdaily4pem${i_pcm_run}.nc
     71mv xoutyearly4pem.nc diags/xoutyearly4pem${i_pcm_run}.nc
    7272cp restartfi.nc starts/restartfi${i_pcm_run}.nc
    7373mv restartfi.nc startfi.nc

  • trunk/LMDZ.COMMON/libf/evolution/deftank/pem_workflow_lib.sh

    r4117 r4152  
    498498    cleanup diags/diagfi .nc $i_resume
    499499    cleanup diags/diagsoil .nc $i_resume
    500     cleanup diags/Xoutdaily4pem .nc $i_resume
    501     cleanup diags/Xoutyearly4pem .nc $i_resume
     500    cleanup diags/xoutdaily4pem .nc $i_resume
     501    cleanup diags/xoutyearly4pem .nc $i_resume
    502502    cleanup logs/run_pcm .log $i_resume
    503503    cleanup starts/restart1D .txt $i_resume
    … …  
    527527        fi
    528528        if [ $i_resume -eq $(($n_pcm_runs_ini - 1)) ]; then
    529             cp diags/Xoutdaily4pem${i_resume}.nc Xoutdaily4pem_Y1.nc
    530             cp diags/Xoutyearly4pem${i_resume}.nc Xoutyearly4pem_Y1.nc
     529            cp diags/xoutdaily4pem${i_resume}.nc xoutdaily4pem_y1.nc
     530            cp diags/xoutyearly4pem${i_resume}.nc xoutyearly4pem_y1.nc
    531531            submit_cycle $1 $n_pcm_runs_ini $i_pcm_run
    532532        elif [ $i_resume -eq $n_pcm_runs_ini ]; then
    533             cp diags/Xoutdaily4pem$(($i_resume - 1)).nc Xoutdaily4pem_Y1.nc
    534             cp diags/Xoutyearly4pem$(($i_resume - 1)).nc Xoutyearly4pem_Y1.nc
    535             cp diags/Xoutdaily4pem${i_resume}.nc Xoutdaily4pem_Y2.nc
    536             cp diags/Xoutyearly4pem${i_resume}.nc Xoutyearly4pem_Y2.nc
     533            cp diags/xoutdaily4pem$(($i_resume - 1)).nc xoutdaily4pem_y1.nc
     534            cp diags/xoutyearly4pem$(($i_resume - 1)).nc xoutyearly4pem_y1.nc
     535            cp diags/xoutdaily4pem${i_resume}.nc xoutdaily4pem_y2.nc
     536            cp diags/xoutyearly4pem${i_resume}.nc xoutyearly4pem_y2.nc
    537537            submit_pem_phase $1 # The next job is a PEM run
    538538        else
    … …  
    554554        cp starts/restartevo$(($i_pem_run - 1)).nc startevo.nc
    555555        if [ $il -eq $(($n_pcm_runs - 1)) ]; then # Second to last PCM run
    556             cp diags/Xoutdaily4pem${i_resume}.nc Xoutdaily4pem_Y1.nc
    557             cp diags/Xoutyearly4pem${i_resume}.nc Xoutyearly4pem_Y1.nc
     556            cp diags/xoutdaily4pem${i_resume}.nc xoutdaily4pem_y1.nc
     557            cp diags/xoutyearly4pem${i_resume}.nc xoutyearly4pem_y1.nc
    558558            submit_cycle $1 $n_pcm_runs $(($il + 1))
    559559        elif [ $il -eq $n_pcm_runs ]; then # Last PCM run so the next job is a PEM run
    560             cp diags/Xoutdaily4pem$(($i_resume - 1)).nc Xoutdaily4pem_Y1.nc
    561             cp diags/Xoutyearly4pem$(($i_resume - 1)).nc Xoutyearly4pem_Y1.nc
    562             cp diags/Xoutdaily4pem${i_resume}.nc Xoutdaily4pem_Y2.nc
    563             cp diags/Xoutyearly4pem${i_resume}.nc Xoutyearly4pem_Y2.nc
     560            cp diags/xoutdaily4pem$(($i_resume - 1)).nc xoutdaily4pem_y1.nc
     561            cp diags/xoutyearly4pem$(($i_resume - 1)).nc xoutyearly4pem_y1.nc
     562            cp diags/xoutdaily4pem${i_resume}.nc xoutdaily4pem_y2.nc
     563            cp diags/xoutyearly4pem${i_resume}.nc xoutyearly4pem_y2.nc
    564564            submit_pem_phase $1
    565565        else
    … …  
    582582    cleanup starts/restart .nc $(($i_pcm_run - 1))
    583583    cleanup starts/restartfi .nc $(($i_pcm_run - 1))
    584     cleanup diags/Xoutdaily4pem .nc $(($i_pcm_run - 1))
    585     cleanup diags/Xoutyearly4pem .nc $(($i_pcm_run - 1))
     584    cleanup diags/xoutdaily4pem .nc $(($i_pcm_run - 1))
     585    cleanup diags/xoutyearly4pem .nc $(($i_pcm_run - 1))
    586586    cleanup diags/diagevo .nc $i_resume
    587587    cleanup diags/diagevo_soil .nc $i_resume

  • trunk/LMDZ.COMMON/libf/evolution/geometry.F90

    r4147 r4152  
    362362!     The longitudes -180 and +180 are duplicated (PCM dynamics).
    363363!     For extensive variables, the values at the poles must be averaged
    364 !     over the number of longitudes for computation on the dynamical grid.
     364!     over the number of longitudes to be computed/outputted on the
     365!     dynamical grid.
    365366!-----------------------------------------------------------------------
    366367

  • trunk/LMDZ.COMMON/libf/evolution/io_netcdf.F90

    r4145 r4152  
    3333character(10), parameter :: startfi_name   = 'startfi.nc'
    3434character(11), parameter :: startevo_name  = 'startevo.nc'
    35 character(19), parameter :: xios_day_name1 = 'Xoutdaily4pem_Y1.nc'  ! XIOS daily output file, second to last PCM year
    36 character(19), parameter :: xios_day_name2 = 'Xoutdaily4pem_Y2.nc'  ! XIOS daily output file, last PCM year
    37 character(20), parameter :: xios_yr_name1  = 'Xoutyearly4pem_Y1.nc' ! XIOS yearly output file, second to last PCM year
    38 character(20), parameter :: xios_yr_name2  = 'Xoutyearly4pem_Y2.nc' ! XIOS yearly output file, last PCM year
     35character(19), parameter :: xios_day_name1 = 'xoutdaily4pem_y1.nc'  ! XIOS daily output file, second to last PCM year
     36character(19), parameter :: xios_day_name2 = 'xoutdaily4pem_y2.nc'  ! XIOS daily output file, last PCM year
     37character(20), parameter :: xios_yr_name1  = 'xoutyearly4pem_y1.nc' ! XIOS yearly output file, second to last PCM year
     38character(20), parameter :: xios_yr_name2  = 'xoutyearly4pem_y2.nc' ! XIOS yearly output file, last PCM year
    3939character(10), parameter :: diagevo_name   = 'diagevo.nc'
    4040

  • trunk/LMDZ.COMMON/libf/evolution/output.F90

    r4147 r4152  
    384384allocate(var_dyn(nlon_loc,nlat))
    385385call vect2dyngrd(longitudes,var_dyn)
    386 call put_var_nc('longitude',var_dyn)
     386call put_var_nc('longitude',var_dyn(:,1))
    387387call vect2dyngrd(latitudes,var_dyn)
    388 call put_var_nc('latitude',var_dyn)
     388call put_var_nc('latitude',var_dyn(1,:))
    389389call put_var_nc('ap',ap)
    390390call put_var_nc('bp',bp)
    391391call put_var_nc('soildepth',mlayer)
    392 call vect2dyngrd(cell_area,var_dyn)
     392call vect2dyngrd(cell_area,var_dyn,extensive = .true.)
    393393call put_var_nc('cell_area',var_dyn)
    394394deallocate(var_dyn)
    … …  
    402402
    403403!=======================================================================
    404 SUBROUTINE write_diagevo(var_name,title,units,var,dimids_in)
     404SUBROUTINE write_diagevo(var_name,title,units,var,dimids,extensive)
    405405!-----------------------------------------------------------------------
    406406! NAME
    … …  
    421421!         - If "diagevo.def" holds a list of variable names, then these
    422422!           are outputted.
    423 !     > If 'var' is an array, then the argument 'dimids_in' is mandatory.
    424 !       'dimids_in' contains the dimension IDs of the variable in order.
     423!     > If 'var' is an array, then the argument 'dimids' is mandatory.
     424!       'dimids' contains the dimension IDs of the variable in order.
    425425!       If the variable holds a dimension 'ngrid', then it must appear
    426426!       in the first position.
     427!     > For extensive variables, the values at the poles must be
     428!       averaged over the number of longitudes to be computed/outputted
     429!       on the dynamical grid.
    427430!-----------------------------------------------------------------------
    428431
    … …  
    442445character(*),               intent(in)           :: var_name, title, units
    443446real(dp),    dimension(..), intent(in)           :: var ! Assumed-rank array
    444 integer(di), dimension(:),  intent(in), optional :: dimids_in
     447integer(di), dimension(:),  intent(in), optional :: dimids
     448logical(k4),                intent(in), optional :: extensive
    445449
    446450! LOCAL VARIABLES
    447451! ---------------
    448 integer(di), dimension(:),       allocatable :: dimids
     452integer(di), dimension(:),       allocatable :: dimids_nc
    449453integer(di)                                  :: idim_ngrid, nlon_loc, ndim, i, j
    450454real(dp),    dimension(:,:),     allocatable :: var_dyn
    … …  
    452456real(dp),    dimension(:,:,:,:), allocatable :: var2d_dyn
    453457integer(di)                                  :: itime ! Current time index to record variables
     458logical(k4)                                  :: extensive_l
    454459
    455460
    … …  
    474479
    475480! Write the variable
     481extensive_l = .false.
     482if (present(extensive)) extensive_l = extensive
    476483call open_nc(diagevo_name,'write') ! Open diagevo_name only it is not done yet
    477484select rank (var) ! Assumed-rank array removes compile-time rank information which is required by the following subroutines using generic interface overload
    … …  
    480487        call put_var_nc(var_name,var,itime)
    481488    rank default ! Arrays
    482         if (.not. present(dimids_in)) call stop_clean(__FILE__,__LINE__,'The argument ''dimids_in'' must be present to output arrays!',1)
    483         ndim = size(dimids_in)
    484         if (any(dimids_in == dim_ngrid)) then
    485             idim_ngrid = findloc(dimids_in,dim_ngrid,dim = 1)
     489        if (.not. present(dimids)) call stop_clean(__FILE__,__LINE__,'The argument ''dimids'' must be present to output arrays!',1)
     490        ndim = size(dimids)
     491        if (any(dimids == dim_ngrid)) then
     492            idim_ngrid = findloc(dimids,dim_ngrid,dim = 1)
    486493            if (idim_ngrid /= 1) call stop_clean(__FILE__,__LINE__,'The dimension ''ngrid'' must be the first one to output the array!',1)
    487             allocate(dimids(ndim + 2))
    488             dimids(1) = dim_nlon
    489             dimids(2) = dim_nlat
    490             dimids(3:ndim + 1) = dimids_in(2:ndim)
    491             dimids(ndim + 2) = dim_time
     494            allocate(dimids_nc(ndim + 2))
     495            dimids_nc(1) = dim_nlon
     496            dimids_nc(2) = dim_nlat
     497            dimids_nc(3:ndim + 1) = dimids(2:ndim)
     498            dimids_nc(ndim + 2) = dim_time
    492499            if (ngrid == 1) then ! 1D
    493500                nlon_loc = 1
    … …  
    495502                nlon_loc = nlon + 1
    496503            end if
    497             call def_var_nc(var_name,title,units,dimids) ! Define the variable only it is not done yet
     504            call def_var_nc(var_name,title,units,dimids_nc) ! Define the variable only it is not done yet
    498505            select rank (var) ! Assumed-rank array removes compile-time rank information which is required by the following subroutines using generic interface overload
    499506                rank(1)
    500507                    allocate(var_dyn(nlon_loc,nlat))
    501                     call vect2dyngrd(var,var_dyn)
     508                    call vect2dyngrd(var,var_dyn,extensive_l)
    502509                    call put_var_nc(var_name,var_dyn,itime)
    503510                    deallocate(var_dyn)
    … …  
    505512                    allocate(var1d_dyn(nlon_loc,nlat,size(var,2)))
    506513                    do i = 1,size(var,2)
    507                         call vect2dyngrd(var(:,i),var1d_dyn(:,:,i))
     514                        call vect2dyngrd(var(:,i),var1d_dyn(:,:,i),extensive_l)
    508515                    end do
    509516                    call put_var_nc(var_name,var1d_dyn,itime)
    … …  
    513520                    do i = 1,size(var,2)
    514521                        do j = 1,size(var,3)
    515                             call vect2dyngrd(var(:,i,j),var2d_dyn(:,:,i,j))
     522                            call vect2dyngrd(var(:,i,j),var2d_dyn(:,:,i,j),extensive_l)
    516523                        end do
    517524                    end do
    … …  
    522529            end select
    523530        else ! No 'ngrid' dimension
    524             allocate(dimids(ndim + 1))
    525             dimids(1:ndim) = dimids_in(:)
    526             dimids(ndim + 1) = dim_time
    527             call def_var_nc(var_name,title,units,dimids) ! Define the variable only it is not done yet
     531            allocate(dimids_nc(ndim + 1))
     532            dimids_nc(1:ndim) = dimids(:)
     533            dimids_nc(ndim + 1) = dim_time
     534            call def_var_nc(var_name,title,units,dimids_nc) ! Define the variable only it is not done yet
    528535            select rank (var) ! Assumed-rank array removes compile-time rank information which is required by the following subroutines using generic interface overload
    529536                rank(1)
    … …  
    539546            end select
    540547        end if
    541         deallocate(dimids)
     548        deallocate(dimids_nc)
    542549end select
    543550call close_nc(diagevo_name)

  • trunk/LMDZ.COMMON/libf/evolution/pem.F90

    r4147 r4152  
    279279idt = 0
    280280do while (n_yr_run < min(nmax_yr_runorb,nmax_yr_run) .and. n_yr_sim < ntot_yr_sim)
    281     call print_msg('**** Iteration of the PEM run [Planetary years]: '//real2str(n_yr_run + dt),LVL_NFO)
     281    call print_msg('**** Iteration of the PEM run [plnt yr]: '//real2str(n_yr_run + dt),LVL_NFO)
    282282    ! Evolve global surface pressure
    283283    call evolve_pressure(d_co2ice,delta_co2_ads,do_sorption,ps_avg_glob_old,ps_avg_glob,ps_avg)

  • trunk/LMDZ.COMMON/libf/evolution/stopping_crit.F90

    r4110 r4152  
    450450        S_atm_2_h2o = S_atm_2_h2o + delta_h2o_ads(i)*cell_area(i)
    451451    else
    452         S_h2o_2_atm = S_h2o_2_atm + delta_h2o_ads(i)*cell_area(i)
     452        S_h2o_2_atm = S_h2o_2_atm - delta_h2o_ads(i)*cell_area(i)
    453453    end if
    454454    if (delta_icetable(i) > 0._dp) then
    455455        S_atm_2_h2o = S_atm_2_h2o + delta_icetable(i)*cell_area(i)
    456456    else
    457         S_h2o_2_atm = S_h2o_2_atm + delta_icetable(i)*cell_area(i)
     457        S_h2o_2_atm = S_h2o_2_atm - delta_icetable(i)*cell_area(i)
    458458    end if
    459459    do islope = 1,nslope

  • trunk/LMDZ.COMMON/libf/evolution/surf_ice.F90

    r4147 r4152  
    378378end if
    379379
    380 h2o_ice = h2o_ice + d_h2oice_new*dt
    381 zshift_surf = d_h2oice_new*dt/rho_h2oice
     380h2o_ice(:,:) = h2o_ice(:,:) + d_h2oice_new(:,:)*dt
     381zshift_surf(:,:) = d_h2oice_new(:,:)*dt/rho_h2oice
    382382
    383383END SUBROUTINE evolve_h2oice
    … …  
    385385
    386386!=======================================================================
    387 SUBROUTINE balance_h2oice_reservoirs(S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
     387SUBROUTINE balance_h2oice_reservoirs(S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_balanced)
    388388!-----------------------------------------------------------------------
    389389! NAME
    … …  
    414414real(dp), dimension(:,:), intent(in)    :: h2o_ice
    415415real(dp), dimension(:,:), intent(inout) :: d_h2oice
    416 real(dp), dimension(:,:), intent(out)   :: d_h2oice_new
     416real(dp), dimension(:,:), intent(out)   :: d_balanced
    417417
    418418! LOCAL VARIABLES
    419419! ---------------
    420420integer(di) :: i, islope
    421 real(qp)    :: S_target, S_target_subl_h2oice, S_target_cond_h2oice, S_ghostice ! Balance variables
     421real(qp)    :: S_corr_subl, S_corr_cond, S_target_subl, S_target_cond, S_ghostice ! Balance variables
    422422real(dp)    :: d_target
    423423
    424424! CODE
    425425! ----
    426 S_target = (S_atm_2_h2o + S_h2o_2_atm)/2._dp
    427 S_target_cond_h2oice = S_atm_2_h2oice + S_target - S_atm_2_h2o
    428 S_target_subl_h2oice = S_h2oice_2_atm + S_target - S_h2o_2_atm
    429 
    430 d_h2oice_new = 0._dp
     426! Compute global targets
     427S_corr_cond = (S_h2o_2_atm - S_atm_2_h2o)/2._qp ! Correction = deviation to the mean
     428S_corr_subl = (S_atm_2_h2o - S_h2o_2_atm)/2._qp ! Correction = deviation to the mean
     429S_target_cond = abs(S_atm_2_h2oice + S_corr_cond)
     430S_target_subl = abs(S_h2oice_2_atm + S_corr_subl)
     431
     432! Compute balanced tendencies with positivity limiter
     433d_balanced(:,:) = 0._dp
    431434S_ghostice = 0._qp
    432435do i = 1,ngrid
    433436    do islope = 1,nslope
    434         if (d_h2oice(i,islope) > 0._dp) then ! Condensing
    435             d_h2oice_new(i,islope) = d_h2oice(i,islope)*real(S_target_cond_h2oice/S_atm_2_h2oice,dp)
    436         else if (d_h2oice(i,islope) < 0._dp) then ! Sublimating
    437             d_target = d_h2oice(i,islope)*real(S_target_subl_h2oice/S_h2oice_2_atm,dp)
    438             if (abs(d_target*dt) <= h2o_ice(i,islope)) then ! Enough ice to sublimate everything
    439                 d_h2oice_new(i,islope) = d_target
    440             else ! Not enough ice to sublimate everything
    441                 ! We sublimate what we can
    442                 d_h2oice_new(i,islope) = h2o_ice(i,islope)/dt
    443                 ! It means the tendency is zero next time
     437        if (d_h2oice(i,islope) > 0._dp) then ! Condensation
     438            d_balanced(i,islope) = d_h2oice(i,islope)*real(S_target_cond/S_atm_2_h2oice,dp)
     439        else if (d_h2oice(i,islope) < 0._dp) then ! Sublimation
     440            d_target = d_h2oice(i,islope)*real(S_target_subl/S_h2oice_2_atm,dp)
     441            if (abs(d_target*dt) <= h2o_ice(i,islope)) then ! Enough ice to sublimate
     442                d_balanced(i,islope) = d_target
     443            else ! Not enough ice to sublimate
     444                ! Sublimate all the available ice
     445                d_balanced(i,islope) = -h2o_ice(i,islope)/dt
     446                ! If fully depleted, zero tendency for next time
    444447                d_h2oice(i,islope) = 0._dp
    445                 ! We compute the amount of H2O ice that we could not make sublimate
    446                 S_ghostice = S_ghostice + abs(d_target*dt) - h2o_ice(i,islope)
     448                ! Compute the amount of H2O ice unable to sublimate
     449                S_ghostice = S_ghostice + real(abs(d_target*dt),qp) - real(h2o_ice(i,islope),qp)
    447450            end if
    448451        end if
    … …  
    450453end do
    451454
    452 ! We need to remove this ice unable to sublimate from places where ice condensed in order to keep balance
    453 where (d_h2oice_new > 0._dp) d_h2oice_new = d_h2oice_new*real((S_target_cond_h2oice - S_ghostice)/S_target_cond_h2oice,dp)
     455! Enforce conservation removing the ghost ice from places where ice condensed
     456where (d_balanced(:,:) > 0._dp) d_balanced(:,:) = d_balanced(:,:)*real((S_target_cond - S_ghostice)/S_target_cond,dp)
    454457
    455458END SUBROUTINE balance_h2oice_reservoirs

  • trunk/LMDZ.COMMON/libf/evolution/tendencies.F90

    r4140 r4152  
    6565
    6666! If the tendency is negative but there is no ice reservoir for the PEM
    67 where (d_ice < 0._dp .and. abs(perice) < minieps) d_ice = 0._dp
     67where (d_ice(:,:) < 0._dp .and. abs(perice(:,:)) < minieps) d_ice(:,:) = 0._dp
    6868
    6969END SUBROUTINE compute_tendice
    … …  
    8989! DEPENDENCIES
    9090! ------------
    91 use geometry,   only: ngrid, nslope, nday
    92 use physics,    only: sigmaB, alpha_clap_co2, beta_clap_co2, Lco2, m_co2, A, B
    93 use orbit,      only: yr_len_sol, sol_len_s
    94 use display,    only: print_msg, LVL_NFO
    95 use utility,    only: real2str
     91use geometry, only: ngrid, nslope, nday
     92use physics,  only: sigmaB, alpha_clap_co2, beta_clap_co2, Lco2, m_co2, A, B
     93use orbit,    only: yr_len_sol, sol_len_s
     94use display,  only: print_msg, LVL_NFO
     95use utility,  only: real2str
    9696
    9797! DECLARATION
Note: See TracChangeset for help on using the changeset viewer.