From 0ca3fa705da1b65892ab32374d15c43a723d5d06 Mon Sep 17 00:00:00 2001 From: Tanya Smirnova Date: Fri, 24 Jan 2025 23:13:09 +0000 Subject: [PATCH] Changes to the code related to fractional sea ice in mpas_atmphys_driver_seaice.F, mpas_atmphys_driver_sfclayer.F, mpas_atmphys_lsm_shared.F: 1. Fixed bug that caused MPAS to crash with frac_seaice=false. It is fixed with adding 'if(config_frac_seaice) then' around the code in seaice driver that uses *_sea variables, which are not allocated. For example, in the line: chs_p(i,j) = xice_p(i,j)*chs_p(i,j) + (1._RKIND-xice_p(i,j))*chs_sea(i,j) 2. We've seen a computational problem with the current use of cpm variable in seaice driver: cpm(i) = xice_p(i,j)*cpm(i) + (1._RKIND-xice_p(i,j))*cpm_sea(i,j) which is fixed by treating cpm as all other variables: cpm_p(i,j) = xice_p(i,j)*cpm_p(i,j) + (1._RKIND-xice_p(i,j))*cpm_sea(i,j) 3. We recomment to replace if(xice_p(i,j).ge.xice_threshold .and. xice_p(i,j).le.1._RKIND) with if(xice_p(i,j).ge.xice_threshold .and. xice_p(i,j).lt.1._RKIND) then to avoid unnecessary computations over vast areas in Arctic covered with sea ice where xice_p(i,j) = 1. --- .../physics/mpas_atmphys_driver_seaice.F | 15 +++++++++++---- .../physics/mpas_atmphys_driver_sfclayer.F | 8 ++++---- .../physics/mpas_atmphys_lsm_shared.F | 2 +- 3 files changed, 16 insertions(+), 9 deletions(-) diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_seaice.F b/src/core_atmosphere/physics/mpas_atmphys_driver_seaice.F index 7894a81429..0b2ab46b37 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_seaice.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_seaice.F @@ -217,6 +217,7 @@ subroutine seaice_from_MPAS(configs,diag_physics,sfc_input,its,ite) qgh_p(i,j) = qgh(i) snoalb_p(i,j) = snoalb(i) br_p(i,j) = br(i) + cpm_p(i,j) = cpm(i) chs_p(i,j) = chs(i) swdown_p(i,j) = gsw(i)/(1._RKIND-sfc_albedo(i)) @@ -252,7 +253,7 @@ subroutine seaice_from_MPAS(configs,diag_physics,sfc_input,its,ite) noahres_p(i,j) = noahres(i) !modify the surface albedo and surface emissivity, and surface temperatures over sea-ice points: - if(xice(i).ge.xice_threshold .and. xice(i).le.1._RKIND) then + if(xice(i).ge.xice_threshold .and. xice(i).lt.1._RKIND) then sfc_albedo_p(i,j) = (sfc_albedo(i) - 0.08_RKIND*(1._RKIND-xice(i))) / xice(i) sfc_emiss_p(i,j) = (sfc_emiss(i) - 0.98_RKIND*(1._RKIND-xice(i))) / xice(i) else @@ -315,10 +316,14 @@ subroutine seaice_to_MPAS(configs,diag_physics,sfc_input,its,ite) !local variables and arrays: integer:: i,j,n +!local pointers: + logical,pointer:: config_frac_seaice !----------------------------------------------------------------------------------------------------------------- !call mpas_log_write('--- enter subroutine seaice_to_MPAS:') + call mpas_pool_get_config(configs,'config_frac_seaice',config_frac_seaice) + call mpas_pool_get_array(diag_physics,'acsnom' ,acsnom ) call mpas_pool_get_array(diag_physics,'acsnow' ,acsnow ) call mpas_pool_get_array(diag_physics,'chs' ,chs ) @@ -351,11 +356,11 @@ subroutine seaice_to_MPAS(configs,diag_physics,sfc_input,its,ite) call mpas_pool_get_array(sfc_input,'xice' ,xice ) !--- weigh local variables needed in the calculation of t2m, th2m, and q2 over seaice points: +if(config_frac_seaice) then do j = jts,jte do i = its,ite - if(xice_p(i,j).ge.xice_threshold .and. xice_p(i,j).le.1._RKIND) then - cpm(i) = xice_p(i,j)*cpm(i) + (1._RKIND-xice_p(i,j))*cpm_sea(i,j) - + if(xice_p(i,j).ge.xice_threshold .and. xice_p(i,j).lt.1._RKIND) then + cpm_p(i,j) = xice_p(i,j)*cpm_p(i,j) + (1._RKIND-xice_p(i,j))*cpm_sea(i,j) chs_p(i,j) = xice_p(i,j)*chs_p(i,j) + (1._RKIND-xice_p(i,j))*chs_sea(i,j) chs2_p(i,j) = xice_p(i,j)*chs2_p(i,j) + (1._RKIND-xice_p(i,j))*chs2_sea(i,j) cqs2_p(i,j) = xice_p(i,j)*cqs2_p(i,j) + (1._RKIND-xice_p(i,j))*cqs2_sea(i,j) @@ -370,6 +375,7 @@ subroutine seaice_to_MPAS(configs,diag_physics,sfc_input,its,ite) endif enddo enddo +endif call sfcdiags( & hfx = hfx_p , qfx = qfx_p , tsk = tsk_p , qsfc = qsfc_p , chs = chs_p , & @@ -409,6 +415,7 @@ subroutine seaice_to_MPAS(configs,diag_physics,sfc_input,its,ite) chs(i) = chs_p(i,j) chs2(i) = chs2_p(i,j) cqs2(i) = cqs2_p(i,j) + cpm(i) = cpm_p(i,j) qsfc(i) = qsfc_p(i,j) qgh(i) = qgh_p(i,j) hfx(i) = hfx_p(i,j) diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F b/src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F index afde4fa523..62aec8c4d2 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F @@ -517,7 +517,7 @@ subroutine sfclayer_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite) v10_sea(i,j) = 0._RKIND !overwrite some local variables for sea-ice cells: - if(xice_p(i,j).ge.xice_threshold .and. xice_p(i,j).le.1._RKIND) then + if(xice_p(i,j).ge.xice_threshold .and. xice_p(i,j).lt.1._RKIND) then xland_sea(i,j) = 2._RKIND mavail_sea(i,j) = 1._RKIND znt_sea(i,j) = 0.0001_RKIND @@ -716,7 +716,7 @@ subroutine sfclayer_to_MPAS(configs,sfc_input,diag_physics,its,ite) call mpas_pool_get_array(sfc_input,'xice',xice) do j = jts,jte do i = its,ite - if(xice(i).ge.xice_threshold .and. xice(i).le.1._RKIND) then + if(xice(i).ge.xice_threshold .and. xice(i).lt.1._RKIND) then br(i) = br_p(i,j)*xice(i) + (1._RKIND-xice(i))*br_sea(i,j) flhc(i) = flhc_p(i,j)*xice(i) + (1._RKIND-xice(i))*flhc_sea(i,j) flqc(i) = flqc_p(i,j)*xice(i) + (1._RKIND-xice(i))*flqc_sea(i,j) @@ -762,7 +762,7 @@ subroutine sfclayer_to_MPAS(configs,sfc_input,diag_physics,its,ite) call mpas_pool_get_array(sfc_input,'xice',xice) do j = jts,jte do i = its,ite - if(xice(i).ge.xice_threshold .and. xice(i).le.1._RKIND) then + if(xice(i).ge.xice_threshold .and. xice(i).lt.1._RKIND) then fh(i) = fh_p(i,j)*xice(i) + (1._RKIND-xice(i))*fh_sea(i,j) fm(i) = fm_p(i,j)*xice(i) + (1._RKIND-xice(i))*fm_sea(i,j) endif @@ -782,7 +782,7 @@ subroutine sfclayer_to_MPAS(configs,sfc_input,diag_physics,its,ite) call mpas_pool_get_array(sfc_input,'xice',xice) do j = jts,jte do i = its,ite - if(xice(i).ge.xice_threshold .and. xice(i).le.1._RKIND) then + if(xice(i).ge.xice_threshold .and. xice(i).lt.1._RKIND) then ch(i) = ch_p(i,j)*xice(i) + (1._RKIND-xice(i))*ch_sea(i,j) endif enddo diff --git a/src/core_atmosphere/physics/mpas_atmphys_lsm_shared.F b/src/core_atmosphere/physics/mpas_atmphys_lsm_shared.F index af5a1a436a..8906f66713 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_lsm_shared.F +++ b/src/core_atmosphere/physics/mpas_atmphys_lsm_shared.F @@ -42,7 +42,7 @@ subroutine correct_tsk_over_seaice(ims,ime,jms,jme,its,ite,jts,jte,xice_thresh,x tsk_sea(i,j) = tsk(i,j) tsk_ice(i,j) = tsk(i,j) - if(xice(i,j).ge.xice_thresh .and. xice(i,j).le.1._RKIND) then + if(xice(i,j).ge.xice_thresh .and. xice(i,j).lt.1._RKIND) then !over sea-ice grid cells, limit sea-surface temperatures to temperatures warmer than 271.4: tsk_sea(i,j) = max(tsk_sea(i,j),271.4_RKIND)