From d6e0bf606cd62896dd8ea4cdd2586be5cb5edc74 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Mon, 27 Jan 2025 09:27:15 -0500 Subject: [PATCH 1/5] v2: Support for MOM6 as of 2024-Nov-27 --- MOM6_GEOSPlug/mom6_cmake/CMakeLists.txt | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/MOM6_GEOSPlug/mom6_cmake/CMakeLists.txt b/MOM6_GEOSPlug/mom6_cmake/CMakeLists.txt index 41ee84e..a3b9b8e 100644 --- a/MOM6_GEOSPlug/mom6_cmake/CMakeLists.txt +++ b/MOM6_GEOSPlug/mom6_cmake/CMakeLists.txt @@ -69,7 +69,9 @@ list( APPEND MOM6_SRCS src/core/MOM_verticalGrid.F90 src/core/MOM_porous_barriers.F90 src/diagnostics/MOM_debugging.F90 + src/diagnostics/MOM_diagnose_MLD.F90 src/diagnostics/MOM_diagnostics.F90 + src/diagnostics/MOM_harmonic_analysis.F90 src/diagnostics/MOM_obsolete_diagnostics.F90 src/diagnostics/MOM_obsolete_params.F90 src/diagnostics/MOM_PointAccel.F90 @@ -176,6 +178,7 @@ list( APPEND MOM6_SRCS src/parameterizations/lateral/MOM_MEKE_types.F90 src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 src/parameterizations/lateral/MOM_self_attr_load.F90 + src/parameterizations/lateral/MOM_streaming_filter.F90 src/parameterizations/lateral/MOM_thickness_diffuse.F90 src/parameterizations/lateral/MOM_spherical_harmonics.F90 src/parameterizations/lateral/MOM_tidal_forcing.F90 @@ -203,7 +206,7 @@ list( APPEND MOM6_SRCS src/parameterizations/vertical/MOM_sponge.F90 src/parameterizations/vertical/MOM_tidal_mixing.F90 src/parameterizations/vertical/MOM_vert_friction.F90 - src/parameterizations/stochastic/MOM_stochastics.F90 + src/parameterizations/stochastic/MOM_stochastics.F90 src/tracer/advection_test_tracer.F90 src/tracer/boundary_impulse_tracer.F90 src/tracer/DOME_tracer.F90 @@ -526,7 +529,7 @@ esma_add_library (${this} SRCS ${SRCS} DEPENDENCIES fms_r8 INCLUDES - $ + $ # choose above set MOM6_infra interface $ $ From 10c6b963d4485e2f799e18fdcf6893615e71b0a5 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Mon, 27 Jan 2025 10:01:14 -0500 Subject: [PATCH 2/5] Update CI and Github Actions --- .circleci/config.yml | 8 ++++---- .github/workflows/enforce-labels.yml | 9 ++++++--- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/.circleci/config.yml b/.circleci/config.yml index 49e3b16..7c7d0d6 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -1,11 +1,11 @@ version: 2.1 -# Anchors to prevent forgetting to update a version -#baselibs_version: &baselibs_version v7.17.0 -#bcs_version: &bcs_version v11.3.0 +# Anchors in case we need to override the defaults from the orb +#baselibs_version: &baselibs_version v7.27.0 +#bcs_version: &bcs_version v11.6.0 orbs: - ci: geos-esm/circleci-tools@2 + ci: geos-esm/circleci-tools@4 workflows: build-test: diff --git a/.github/workflows/enforce-labels.yml b/.github/workflows/enforce-labels.yml index 6e1720e..86f4bb4 100644 --- a/.github/workflows/enforce-labels.yml +++ b/.github/workflows/enforce-labels.yml @@ -8,18 +8,20 @@ jobs: require-label: runs-on: ubuntu-latest steps: - - uses: mheap/github-action-required-labels@v2 + - uses: mheap/github-action-required-labels@v5 env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} with: mode: minimum count: 1 - labels: "0 diff,0 diff trivial,Non 0-diff,0 diff structural,0-diff trivial,Not 0-diff,0-diff,automatic,0-diff uncoupled" + labels: "0 diff,0 diff trivial,Non 0-diff,0 diff structural,0-diff trivial,Not 0-diff,0-diff,automatic,0-diff uncoupled,github_actions" add_comment: true + message: "This PR is being prevented from merging because you have not added one of our required labels: {{ provided }}. Please add one so that the PR can be merged." + blocking-label: runs-on: ubuntu-latest steps: - - uses: mheap/github-action-required-labels@v2 + - uses: mheap/github-action-required-labels@v5 env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} with: @@ -27,3 +29,4 @@ jobs: count: 0 labels: "Contingent - DNA,Needs Lead Approval,Contingent -- Do Not Approve" add_comment: true + message: "This PR is being prevented from merging because you have added one of our blocking labels: {{ provided }}. You'll need to remove it before this PR can be merged." From 03dbeda06fa3627b2486aeb073f20defdaf5cc02 Mon Sep 17 00:00:00 2001 From: bzhao Date: Fri, 14 Feb 2025 10:21:58 -0500 Subject: [PATCH 3/5] fixed a bug computing FRZMLT --- MOM6_GEOSPlug/MOM6_GEOSPlug.F90 | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/MOM6_GEOSPlug/MOM6_GEOSPlug.F90 b/MOM6_GEOSPlug/MOM6_GEOSPlug.F90 index 58f07b7..b013761 100644 --- a/MOM6_GEOSPlug/MOM6_GEOSPlug.F90 +++ b/MOM6_GEOSPlug/MOM6_GEOSPlug.F90 @@ -830,10 +830,15 @@ subroutine Run ( GC, IMPORT, EXPORT, CLOCK, RC ) endif where(MOM_2D_MASK(:,:)>0.0) - FRZMLT = FRAZIL + MELT_POT + FRZMLT = MELT_POT elsewhere FRZMLT = 0.0 end where + + where(MOM_2D_MASK(:,:)>0.0 .and. FRAZIL(:,:) > 0.0) + FRZMLT = FRAZIL + endwhere + end if ! freezing temperature (deg C) From 0c5411f9df8f40d59e11c8ea31320ce4a97ef070 Mon Sep 17 00:00:00 2001 From: Michael Mehari Date: Thu, 6 Mar 2025 15:58:45 -0500 Subject: [PATCH 4/5] Added export for MOM6 T, S and DH fields to be imported by NOBM --- MOM6_GEOSPlug/MOM6_GEOSPlug.F90 | 90 +++++++++++++++++++---- MOM6_GEOSPlug/MOM6_GEOSPlug_StateSpecs.rc | 8 +- 2 files changed, 78 insertions(+), 20 deletions(-) diff --git a/MOM6_GEOSPlug/MOM6_GEOSPlug.F90 b/MOM6_GEOSPlug/MOM6_GEOSPlug.F90 index b013761..783c40e 100644 --- a/MOM6_GEOSPlug/MOM6_GEOSPlug.F90 +++ b/MOM6_GEOSPlug/MOM6_GEOSPlug.F90 @@ -44,16 +44,19 @@ module MOM6_GEOSPlugMod use MOM_grid, only: ocean_grid_type - use ocean_model_mod, only: get_ocean_grid, & - ocean_model_init, & - ocean_model_init_sfc, & - update_ocean_model, & - ocean_model_end, & - ocean_model_restart, & - ocean_model_data_get, & - ocean_public_type, & - ocean_state_type, & - ocean_model_get_UV_surf + use ocean_model_mod, only: get_ocean_grid, & + ocean_model_init, & + ocean_model_init_sfc, & + update_ocean_model, & + ocean_model_end, & + ocean_model_restart, & + ocean_model_data_get, & + ocean_public_type, & + ocean_state_type, & + ocean_model_get_UV_surf, & + ocean_model_get_thickness, & + ocean_model_get_prog_tracer, & + ocean_model_get_prog_tracer_index use MOM_surface_forcing_gfdl, only: ice_ocean_boundary_type @@ -182,7 +185,7 @@ subroutine Initialize ( GC, IMPORT, EXPORT, CLOCK, RC ) integer :: Comm integer :: isc,iec,jsc,jec integer :: isd,ied,jsd,jed - integer :: IM, JM + integer :: IM, JM, LM integer :: g_isc,g_iec,g_jsc,g_jec integer :: g_isd,g_ied,g_jsd,g_jed integer :: YEAR,MONTH,DAY,HR,MN,SC @@ -216,11 +219,14 @@ subroutine Initialize ( GC, IMPORT, EXPORT, CLOCK, RC ) REAL_, pointer :: MOM_2D_MASK(:,:) => null() REAL_, pointer :: SLV(:,:) => null() - real, allocatable :: Tmp2(:,:) + real, allocatable :: Tmp2(:,:), Tmp3(:,:,:) + + REAL_, pointer, dimension(:, :, :) :: DH, TL, SL integer :: DT_OCEAN character(len=7) :: wind_stagger ! 'AGRID' or 'BGRID' or 'CGRID' integer ::iwind_stagger ! AGRID or BGRID or CGRID + integer :: i __Iam__('Initialize') @@ -353,6 +359,9 @@ subroutine Initialize ( GC, IMPORT, EXPORT, CLOCK, RC ) ASSERT_(counts(1)==IM) ASSERT_(counts(2)==JM) + call get_ocean_grid(Ocean_state, Ocean_grid) + LM=Ocean_grid%ke + ! Check run time surface current stagger option set in MOM_input ! to make sure it matches what is expected here !--------------------------------------------------------------- @@ -443,6 +452,24 @@ subroutine Initialize ( GC, IMPORT, EXPORT, CLOCK, RC ) call MAPL_GetPointer(EXPORT, SW, 'SW' , alloc=.true., _RC) call MAPL_GetPointer(EXPORT, AREA, 'AREA', alloc=.true., _RC) call MAPL_GetPointer(EXPORT, SLV, 'SLV', alloc=.true., _RC) + call MAPL_GetPointer(EXPORT, DH, 'DH' , alloc=.true., _RC) + call MAPL_GetPointer(EXPORT, TL, 'T', alloc=.true., _RC) + call MAPL_GetPointer(EXPORT, SL, 'S', alloc=.true., _RC) + +! Get the 3-D MOM data +!--------------------- + allocate(Tmp3(IM,JM,LM), __STAT__) + + call ocean_model_get_thickness(Ocean_State, tmp3, isc, jsc) + DH = real(tmp3, kind=GeosKind) + + call ocean_model_get_prog_tracer_index(Ocean_State,i,'temp') + call ocean_model_get_prog_tracer(Ocean_State,i, tmp3, isc, jsc) + TL = real(tmp3,kind=GeosKind) + MAPL_TICE + + call ocean_model_get_prog_tracer_index(Ocean_State,i,'salt') + call ocean_model_get_prog_tracer(Ocean_State, i, tmp3, isc, jsc) + SL = real(tmp3,kind=GeosKind) ! Get the 2-D MOM data !--------------------- @@ -471,6 +498,7 @@ subroutine Initialize ( GC, IMPORT, EXPORT, CLOCK, RC ) end if deallocate(Tmp2, __STAT__) + deallocate(Tmp3, __STAT__) RETURN_(ESMF_SUCCESS) end subroutine Initialize @@ -511,6 +539,7 @@ subroutine Run ( GC, IMPORT, EXPORT, CLOCK, RC ) type(ice_ocean_boundary_type), pointer :: Boundary => null() type(ocean_public_type), pointer :: Ocean => null() type(ocean_state_type), pointer :: Ocean_State => null() + type(ocean_grid_type), pointer :: Ocean_grid => null() type(MOM_MAPL_Type), pointer :: MOM_MAPL_internal_state => null() type(MOM_MAPLWrap_Type) :: wrap @@ -531,6 +560,9 @@ subroutine Run ( GC, IMPORT, EXPORT, CLOCK, RC ) REAL_, pointer :: MOM_2D_MASK (:,:) => null() REAL_, pointer :: AREA (:,:) => null() REAL_, pointer :: T_Freeze (:,:) => null() + REAL_, pointer :: DH (:,:,:) => null() + REAL_, pointer :: TL (:,:,:) => null() + REAL_, pointer :: SL (:,:,:) => null() ! Imports REAL_, pointer :: TAUX(:,:) => null() @@ -555,11 +587,11 @@ subroutine Run ( GC, IMPORT, EXPORT, CLOCK, RC ) REAL_, pointer :: TAUYBOT(:,:) => null() ! Temporaries - real, allocatable :: U (:,:), V(:,:) + real, allocatable :: U(:,:), V(:,:), H(:,:,:) real, allocatable :: cos_rot(:,:) real, allocatable :: sin_rot(:,:) - integer :: IM, JM + integer :: IM, JM, LM integer :: steady_state_ocean = 0 ! SA: Per Atanas T, "name" of this var is misleading ! We run ocean model only when it = 0 @@ -575,6 +607,7 @@ subroutine Run ( GC, IMPORT, EXPORT, CLOCK, RC ) real :: pice_scaling = 1.0 real :: dTf_dS integer :: DT_OCEAN + integer :: tracer_index REAL_, pointer, dimension(:,:) :: LATS => null() REAL_, pointer, dimension(:,:) :: LONS => null() @@ -626,11 +659,15 @@ subroutine Run ( GC, IMPORT, EXPORT, CLOCK, RC ) IM=iec-isc+1 JM=jec-jsc+1 + call get_ocean_grid (Ocean_state, Ocean_grid) + LM=Ocean_grid%ke + ! Temporaries with MOM default reals !----------------------------------- allocate(U(IM,JM ), __STAT__) allocate(V(IM,JM ), __STAT__) + allocate(H(IM,JM,LM), __STAT__) allocate(cos_rot(IM,JM), __STAT__) allocate(sin_rot(IM,JM), __STAT__) @@ -662,6 +699,10 @@ subroutine Run ( GC, IMPORT, EXPORT, CLOCK, RC ) ! Get EXPORT pointers !-------------------- + call MAPL_GetPointer(EXPORT, DH, 'DH' , _RC) + call MAPL_GetPointer(EXPORT, TL, 'T' , _RC) + call MAPL_GetPointer(EXPORT, SL, 'S' , _RC) + call MAPL_GetPointer(EXPORT, UW, 'UW' , _RC) call MAPL_GetPointer(EXPORT, VW, 'VW' , _RC) call MAPL_GetPointer(EXPORT, UWB, 'UWB' , _RC) @@ -765,7 +806,23 @@ subroutine Run ( GC, IMPORT, EXPORT, CLOCK, RC ) ! Get required Exports at GEOS precision !--------------------------------------- -! mask +!MOM 3D fields +! Thickness + call ocean_model_get_thickness(Ocean_State, H, isc, jsc) + DH = real(H, kind=GeosKind) + +! Temp + call ocean_model_get_prog_tracer_index(Ocean_State,tracer_index,'temp') + call ocean_model_get_prog_tracer(Ocean_State,tracer_index, H, isc, jsc) + TL = real(H,kind=GeosKind) + MAPL_TICE + +! Salt + call ocean_model_get_prog_tracer_index(Ocean_State,tracer_index,'salt') + call ocean_model_get_prog_tracer(Ocean_State, tracer_index, H, isc, jsc) + SL = real(H,kind=GeosKind) + +!MOM 2D fields +! 2D mask U = 0.0 call ocean_model_data_get(Ocean_State, Ocean, 'mask', U, isc, jsc) MOM_2D_MASK = real(U, kind=GeosKind) @@ -910,11 +967,12 @@ subroutine Run ( GC, IMPORT, EXPORT, CLOCK, RC ) ! Optional Exports at GEOS precision !----------------------------------- ! none -! 3d exports with MOM6, such as depths, T, S, U, V, etc +! 3d exports with MOM6, other than depth, T, S (needed by NOBM), such as U, V, etc ! will not be exported. If needed, write them directly from MOM6 deallocate(U, V, __STAT__) deallocate(cos_rot,sin_rot, __STAT__) + deallocate(H, __STAT__) call MAPL_TimerOff(MAPL,"RUN" ) call MAPL_TimerOff(MAPL,"TOTAL" ) diff --git a/MOM6_GEOSPlug/MOM6_GEOSPlug_StateSpecs.rc b/MOM6_GEOSPlug/MOM6_GEOSPlug_StateSpecs.rc index 30e5fbe..41e7a62 100644 --- a/MOM6_GEOSPlug/MOM6_GEOSPlug_StateSpecs.rc +++ b/MOM6_GEOSPlug/MOM6_GEOSPlug_StateSpecs.rc @@ -49,10 +49,10 @@ SLV | m | xy | N | sea_level_with_ice_loading_and_inv FRAZIL | W m-2 | xy | N | heating_from_frazil_formation MELT_POT | W m-2 | xy | N | heat_that_can_be_used_to_melt_sea_ice FRZMLT | W m-2 | xy | N | freeze_melt_potential -T_Freeze | degC | xy | N | freezing_temperature_calculated_using_salinity -DH | m OR kg m-2 | xyz | C | DUMMY_EXPORT_layer_thickness -T | K | xyz | C | DUMMY_EXPORT_potential_temperature -S | PSU | xyz | C | DUMMY_EXPORT_salinity +T_Freeze | degC | xy | N | freezing_temperature_calculated_using_salinity +DH | m OR kg m-2 | xyz | C | layer_thickness +T | K | xyz | C | potential_temperature +S | PSU | xyz | C | salinity category: INTERNAL #---------------------------------------------------------------------------------------- From dfe0ac30f795bec93e0d74335f7c3f1b586ba23a Mon Sep 17 00:00:00 2001 From: afahadabdullah Date: Fri, 11 Apr 2025 14:22:14 -0700 Subject: [PATCH 5/5] mit C180 LLC270 config --- .../c180_llc270_01/code/CPP_EEOPTIONS.h | 156 ++++ .../configs/c180_llc270_01/code/CPP_OPTIONS.h | 195 ++++ .../c180_llc270_01/code/DIAGNOSTICS_SIZE.h | 28 + .../c180_llc270_01/code/GGL90_OPTIONS.h | 31 + .../c180_llc270_01/code/GMREDI_OPTIONS.h | 69 ++ .../c180_llc270_01/code/MOM_COMMON_OPTIONS.h | 25 + .../c180_llc270_01/code/SEAICE_LAYERS.h | 86 ++ .../c180_llc270_01/code/SEAICE_OPTIONS.h | 264 ++++++ .../configs/c180_llc270_01/code/SEAICE_SIZE.h | 54 ++ .../configs/c180_llc270_01/code/SIZE.h | 65 ++ .../c180_llc270_01/code/mom_calc_visc.F | 794 +++++++++++++++++ .../configs/c180_llc270_01/code/packages.conf | 7 + .../c180_llc270_01/code/seaice_advdiff.F | 771 ++++++++++++++++ .../code/seaice_diag_init_add.h | 120 +++ .../code/seaice_diagnostics_init.F | 838 ++++++++++++++++++ .../c180_llc270_01/code/seaice_model.F | 390 ++++++++ .../c180_llc270_01/code/seaice_save4gmao.F | 261 ++++++ .../configs/c180_llc270_01/input/data | 132 +++ .../c180_llc270_01/input/data.diagnostics | 127 +++ .../configs/c180_llc270_01/input/data.exch2 | 318 +++++++ .../configs/c180_llc270_01/input/data.ggl90 | 15 + .../configs/c180_llc270_01/input/data.gmredi | 38 + .../configs/c180_llc270_01/input/data.pkg | 8 + .../c180_llc270_01/input/data.salt_plume | 6 + .../configs/c180_llc270_01/input/data.seaice | 50 ++ .../configs/c180_llc270_01/input/eedata | 13 + 26 files changed, 4861 insertions(+) create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/CPP_EEOPTIONS.h create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/CPP_OPTIONS.h create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/DIAGNOSTICS_SIZE.h create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/GGL90_OPTIONS.h create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/GMREDI_OPTIONS.h create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/MOM_COMMON_OPTIONS.h create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/SEAICE_LAYERS.h create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/SEAICE_OPTIONS.h create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/SEAICE_SIZE.h create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/SIZE.h create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/mom_calc_visc.F create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/packages.conf create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_advdiff.F create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_diag_init_add.h create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_diagnostics_init.F create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_model.F create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_save4gmao.F create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.diagnostics create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.exch2 create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.ggl90 create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.gmredi create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.pkg create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.salt_plume create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.seaice create mode 100755 MIT_GEOS5PlugMod/configs/c180_llc270_01/input/eedata diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/CPP_EEOPTIONS.h b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/CPP_EEOPTIONS.h new file mode 100755 index 0000000..883fe0d --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/CPP_EEOPTIONS.h @@ -0,0 +1,156 @@ +#ifndef _CPP_EEOPTIONS_H_ +#define _CPP_EEOPTIONS_H_ + +CBOP +C !ROUTINE: CPP_EEOPTIONS.h +C !INTERFACE: +C include "CPP_EEOPTIONS.h" +C +C !DESCRIPTION: +C *==========================================================* +C | CPP\_EEOPTIONS.h | +C *==========================================================* +C | C preprocessor "execution environment" supporting | +C | flags. Use this file to set flags controlling the | +C | execution environment in which a model runs - as opposed | +C | to the dynamical problem the model solves. | +C | Note: Many options are implemented with both compile time| +C | and run-time switches. This allows options to be | +C | removed altogether, made optional at run-time or | +C | to be permanently enabled. This convention helps | +C | with the data-dependence analysis performed by the | +C | adjoint model compiler. This data dependency | +C | analysis can be upset by runtime switches that it | +C | is unable to recoginise as being fixed for the | +C | duration of an integration. | +C | A reasonable way to use these flags is to | +C | set all options as selectable at runtime but then | +C | once an experimental configuration has been | +C | identified, rebuild the code with the appropriate | +C | options set at compile time. | +C *==========================================================* +CEOP + +C In general the following convention applies: +C ALLOW - indicates an feature will be included but it may +C CAN have a run-time flag to allow it to be switched +C on and off. +C If ALLOW or CAN directives are "undef'd" this generally +C means that the feature will not be available i.e. it +C will not be included in the compiled code and so no +C run-time option to use the feature will be available. +C +C ALWAYS - indicates the choice will be fixed at compile time +C so no run-time option will be present + +C=== Macro related options === +C-- Control storage of floating point operands +C On many systems it improves performance only to use +C 8-byte precision for time stepped variables. +C Constant in time terms ( geometric factors etc.. ) +C can use 4-byte precision, reducing memory utilisation and +C boosting performance because of a smaller working set size. +C However, on vector CRAY systems this degrades performance. +C Enable to switch REAL4_IS_SLOW from genmake2 (with LET_RS_BE_REAL4): +#ifdef LET_RS_BE_REAL4 +#undef REAL4_IS_SLOW +#else /* LET_RS_BE_REAL4 */ +#define REAL4_IS_SLOW +#endif /* LET_RS_BE_REAL4 */ + +C-- Control use of "double" precision constants. +C Use D0 where it means REAL*8 but not where it means REAL*16 +#define D0 d0 + +C=== IO related options === +C-- Flag used to indicate whether Fortran formatted write +C and read are threadsafe. On SGI the routines can be thread +C safe, on Sun it is not possible - if you are unsure then +C undef this option. +#undef FMTFTN_IO_THREAD_SAFE + +C-- Flag used to indicate whether Binary write to Local file (i.e., +C a different file for each tile) and read are thread-safe. +#undef LOCBIN_IO_THREAD_SAFE + +C-- Flag to turn off the writing of error message to ioUnit zero +#undef DISABLE_WRITE_TO_UNIT_ZERO + +C-- Alternative formulation of BYTESWAP, faster than +C compiler flag -byteswapio on the Altix. +#undef FAST_BYTESWAP + +C-- Flag to turn on old default of opening scratch files with the +C STATUS='SCRATCH' option. This method, while perfectly FORTRAN-standard, +C caused filename conflicts on some multi-node/multi-processor platforms +C in the past and has been replace by something (hopefully) more robust. +#undef USE_FORTRAN_SCRATCH_FILES + +C-- Flag defined for eeboot_minimal.F, eeset_parms.F and open_copy_data_file.F +C to write STDOUT, STDERR and scratch files from process 0 only. +C WARNING: to use only when absolutely confident that the setup is working +C since any message (error/warning/print) from any proc <> 0 will be lost. +#undef SINGLE_DISK_IO + +C=== MPI, EXCH and GLOBAL_SUM related options === +C-- Flag turns off MPI_SEND ready_to_receive polling in the +C gather_* subroutines to speed up integrations. +#undef DISABLE_MPI_READY_TO_RECEIVE + +C-- Control MPI based parallel processing +CXXX We no longer select the use of MPI via this file (CPP_EEOPTIONS.h) +CXXX To use MPI, use an appropriate genmake2 options file or use +CXXX genmake2 -mpi . +CXXX #undef ALLOW_USE_MPI + +C-- Control use of communication that might overlap computation. +C Under MPI selects/deselects "non-blocking" sends and receives. +#undef ALLOW_ASYNC_COMMUNICATION +#undef ALWAYS_USE_ASYNC_COMMUNICATION +C-- Control use of communication that is atomic to computation. +C Under MPI selects/deselects "blocking" sends and receives. +#define ALLOW_SYNC_COMMUNICATION +#undef ALWAYS_USE_SYNC_COMMUNICATION + +C-- Control XY periodicity in processor to grid mappings +C Note: Model code does not need to know whether a domain is +C periodic because it has overlap regions for every box. +C Model assume that these values have been +C filled in some way. +#undef ALWAYS_PREVENT_X_PERIODICITY +#undef ALWAYS_PREVENT_Y_PERIODICITY +#define CAN_PREVENT_X_PERIODICITY +#define CAN_PREVENT_Y_PERIODICITY + +C-- disconnect tiles (no exchange between tiles, just fill-in edges +C assuming locally periodic subdomain) +#undef DISCONNECTED_TILES + +C-- Always cumulate tile local-sum in the same order by applying MPI allreduce +C to array of tiles ; can get slower with large number of tiles (big set-up) +#define GLOBAL_SUM_ORDER_TILES + +C-- Alternative way of doing global sum without MPI allreduce call +C but instead, explicit MPI send & recv calls. Expected to be slower. +#undef GLOBAL_SUM_SEND_RECV + +C-- Alternative way of doing global sum on a single CPU +C to eliminate tiling-dependent roundoff errors. Note: This is slow. +#undef CG2D_SINGLECPU_SUM + +C=== Other options (to add/remove pieces of code) === +C-- Flag to turn on checking for errors from all threads and procs +C (calling S/R STOP_IF_ERROR) before stopping. +#define USE_ERROR_STOP + +C-- Control use of communication with other component: +C allow to import and export from/to Coupler interface. +#undef COMPONENT_MODULE + +C-- Activate some pieces of code for coupling to GEOS AGCM +#define HACK_FOR_GMAO_CPL + +C=== And define Macros === +#include "CPP_EEMACROS.h" + +#endif /* _CPP_EEOPTIONS_H_ */ diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/CPP_OPTIONS.h b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/CPP_OPTIONS.h new file mode 100755 index 0000000..2a9d70d --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/CPP_OPTIONS.h @@ -0,0 +1,195 @@ +#ifndef CPP_OPTIONS_H +#define CPP_OPTIONS_H + +CBOP +C !ROUTINE: CPP_OPTIONS.h +C !INTERFACE: +C #include "CPP_OPTIONS.h" + +C !DESCRIPTION: +C *==================================================================* +C | main CPP options file for the model: +C | Control which optional features to compile in model/src code. +C *==================================================================* +CEOP + +C CPP flags controlling particular source code features + +C-- Forcing code options: + +C o Shortwave heating as extra term in external_forcing.F +C Note: this should be a run-time option +#define SHORTWAVE_HEATING + +C o Include/exclude Geothermal Heat Flux at the bottom of the ocean +#define ALLOW_GEOTHERMAL_FLUX + +C o Allow to account for heating due to friction (and momentum dissipation) +#undef ALLOW_FRICTION_HEATING + +C o Allow mass source or sink of Fluid in the interior +C (3-D generalisation of oceanic real-fresh water flux) +#undef ALLOW_ADDFLUID + +C o Include pressure loading code +#define ATMOSPHERIC_LOADING + +C o Include/exclude balancing surface forcing fluxes code +#undef ALLOW_BALANCE_FLUXES + +C o Include/exclude balancing surface forcing relaxation code +#undef ALLOW_BALANCE_RELAX + +C o Include/exclude checking for negative salinity +#undef CHECK_SALINITY_FOR_NEGATIVE_VALUES + +C-- Options to discard parts of the main code: + +C o Exclude/allow external forcing-fields load +C this allows to read & do simple linear time interpolation of oceanic +C forcing fields, if no specific pkg (e.g., EXF) is used to compute them. +#undef EXCLUDE_FFIELDS_LOAD +C If defined, use same method (with pkg/autodiff compiled or not) for checking +C when to load new reccord ; by default, use simpler method with pkg/autodiff. +#undef STORE_LOADEDREC_TEST + +C o Include/exclude phi_hyd calculation code +#define INCLUDE_PHIHYD_CALCULATION_CODE + +C o Include/exclude sound speed calculation code +C o (Note that this is a diagnostic from Del Grasso algorithm, not derived +C from EOS) +#undef INCLUDE_SOUNDSPEED_CALC_CODE + +C-- Vertical mixing code options: + +C o Include/exclude calling S/R CONVECTIVE_ADJUSTMENT +#define INCLUDE_CONVECT_CALL + +C o Include/exclude calling S/R CONVECTIVE_ADJUSTMENT_INI, turned off by +C default because it is an unpopular historical left-over +#undef INCLUDE_CONVECT_INI_CALL + +C o Include/exclude call to S/R CALC_DIFFUSIVITY +#define INCLUDE_CALC_DIFFUSIVITY_CALL + +C o Allow full 3D specification of vertical diffusivity +#define ALLOW_3D_DIFFKR + +C o Allow latitudinally varying BryanLewis79 vertical diffusivity +#undef ALLOW_BL79_LAT_VARY + +C o Exclude/allow partial-cell effect (physical or enhanced) in vertical mixing +C this allows to account for partial-cell in vertical viscosity and diffusion, +C either from grid-spacing reduction effect or as artificially enhanced mixing +C near surface & bottom for too thin grid-cell +#undef EXCLUDE_PCELL_MIX_CODE + +C o Exclude/allow to use isotropic 3-D Smagorinsky viscosity as diffusivity +C for tracers (after scaling by constant Prandtl number) +#undef ALLOW_SMAG_3D_DIFFUSIVITY + +C-- Time-stepping code options: + +C o Include/exclude combined Surf.Pressure and Drag Implicit solver code +#undef ALLOW_SOLVE4_PS_AND_DRAG + +C o Include/exclude Implicit vertical advection code +#define INCLUDE_IMPLVERTADV_CODE + +C o Include/exclude AdamsBashforth-3rd-Order code +#undef ALLOW_ADAMSBASHFORTH_3 + +C o Include/exclude Quasi-Hydrostatic Stagger Time-step AdamsBashforth code +#undef ALLOW_QHYD_STAGGER_TS + +C-- Model formulation options: + +C o Allow/exclude "Exact Convervation" of fluid in Free-Surface formulation +C that ensures that d/dt(eta) is exactly equal to - Div.Transport +#define EXACT_CONSERV + +C o Allow the use of Non-Linear Free-Surface formulation +C this implies that grid-cell thickness (hFactors) varies with time +#define NONLIN_FRSURF +C o Disable code for rStar coordinate and/or code for Sigma coordinate +c#define DISABLE_RSTAR_CODE +#define DISABLE_SIGMA_CODE + +C o Include/exclude nonHydrostatic code +#undef ALLOW_NONHYDROSTATIC + +C o Include/exclude GM-like eddy stress in momentum code +#undef ALLOW_EDDYPSI + +C-- Algorithm options: + +C o Include/exclude code for Non Self-Adjoint (NSA) conjugate-gradient solver +#undef ALLOW_CG2D_NSA + +C o Include/exclude code for single reduction Conjugate-Gradient solver +#undef ALLOW_SRCG + +C o Choices for implicit solver routines solve_*diagonal.F +C The following has low memory footprint, but not suitable for AD +#undef SOLVE_DIAGONAL_LOWMEMORY +C The following one suitable for AD but does not vectorize +#define SOLVE_DIAGONAL_KINNER + +C Implementation alternative (might be faster on some platforms ?) +#undef USE_MASK_AND_NO_IF + +C-- Retired code options: + +C o ALLOW isotropic scaling of harmonic and bi-harmonic terms when +C using an locally isotropic spherical grid with (dlambda) x (dphi*cos(phi)) +C *only for use on a lat-lon grid* +C Setting this flag here affects both momentum and tracer equation unless +C it is set/unset again in other header fields (e.g., GAD_OPTIONS.h). +C The definition of the flag is commented to avoid interference with +C such other header files. +C The preferred method is specifying a value for viscAhGrid or viscA4Grid +C in data which is then automatically scaled by the grid size; +C the old method of specifying viscAh/viscA4 and this flag is provided +C for completeness only (and for use with the adjoint). +c#define ISOTROPIC_COS_SCALING + +C o This flag selects the form of COSINE(lat) scaling of bi-harmonic term. +C *only for use on a lat-lon grid* +C Has no effect if ISOTROPIC_COS_SCALING is undefined. +C Has no effect on vector invariant momentum equations. +C Setting this flag here affects both momentum and tracer equation unless +C it is set/unset again in other header fields (e.g., GAD_OPTIONS.h). +C The definition of the flag is commented to avoid interference with +C such other header files. +c#define COSINEMETH_III + +C o Use "OLD" UV discretisation near boundaries (*not* recommended) +C Note - only works with pkg/mom_fluxform and "no_slip_sides=.FALSE." +C because the old code did not have no-slip BCs +#undef OLD_ADV_BCS + +C o Use LONG.bin, LATG.bin, etc., initialization for ini_curviliear_grid.F +C Default is to use "new" grid files (OLD_GRID_IO undef) but OLD_GRID_IO +C is still useful with, e.g., single-domain curvilinear configurations. +#undef OLD_GRID_IO + +C o Use old EXTERNAL_FORCING_U,V,T,S subroutines (for backward compatibility) +#undef USE_OLD_EXTERNAL_FORCING + +C-- Other option files: + +C o Execution environment support options +#include "CPP_EEOPTIONS.h" + +C o Include/exclude single header file containing multiple packages options +C (AUTODIFF, COST, CTRL, ECCO, EXF ...) instead of the standard way where +C each of the above pkg get its own options from its specific option file. +C Although this method, inherited from ECCO setup, has been traditionally +C used for all adjoint built, work is in progress to allow to use the +C standard method also for adjoint built. +c#ifdef PACKAGES_CONFIG_H +c# include "ECCO_CPPOPTIONS.h" +c#endif + +#endif /* CPP_OPTIONS_H */ diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/DIAGNOSTICS_SIZE.h b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/DIAGNOSTICS_SIZE.h new file mode 100755 index 0000000..5ebff97 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/DIAGNOSTICS_SIZE.h @@ -0,0 +1,28 @@ +C Diagnostics Array Dimension +C --------------------------- +C ndiagMax :: maximum total number of available diagnostics +C numlists :: maximum number of diagnostics list (in data.diagnostics) +C numperlist :: maximum number of active diagnostics per list (data.diagnostics) +C numLevels :: maximum number of levels to write (data.diagnostics) +C numDiags :: maximum size of the storage array for active 2D/3D diagnostics +C nRegions :: maximum number of regions (statistics-diagnostics) +C sizRegMsk :: maximum size of the regional-mask (statistics-diagnostics) +C nStats :: maximum number of statistics (e.g.: aver,min,max ...) +C diagSt_size:: maximum size of the storage array for statistics-diagnostics +C Note : may need to increase "numDiags" when using several 2D/3D diagnostics, +C and "diagSt_size" (statistics-diags) since values here are deliberately small. + INTEGER ndiagMax + INTEGER numlists, numperlist, numLevels + INTEGER numDiags + INTEGER nRegions, sizRegMsk, nStats + INTEGER diagSt_size + PARAMETER( ndiagMax = 700 ) + PARAMETER( numlists = 300, numperlist = 30, numLevels=5*Nr ) + PARAMETER( numDiags = 3000 ) + PARAMETER( nRegions = 20 , sizRegMsk = 1 , nStats = 4 ) + PARAMETER( diagSt_size = 50*Nr ) + + +CEH3 ;;; Local Variables: *** +CEH3 ;;; mode:fortran *** +CEH3 ;;; End: *** diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/GGL90_OPTIONS.h b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/GGL90_OPTIONS.h new file mode 100755 index 0000000..64ecd78 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/GGL90_OPTIONS.h @@ -0,0 +1,31 @@ +C *=============================================================* +C | GGL90_OPTIONS.h +C | o CPP options file for GGL90 package. +C *=============================================================* +C | Use this file for selecting options within the GGL90 +C | package. +C *=============================================================* + +#ifndef GGL90_OPTIONS_H +#define GGL90_OPTIONS_H +#include "PACKAGES_CONFIG.h" +#include "CPP_OPTIONS.h" + +#ifdef ALLOW_GGL90 +C Package-specific Options & Macros go here + +C Enable horizontal diffusion of TKE. +#undef ALLOW_GGL90_HORIZDIFF + +C Use horizontal averaging for viscosity and diffusivity as +C originally implemented in OPA. +#define ALLOW_GGL90_SMOOTH + +C allow IDEMIX model +#undef ALLOW_GGL90_IDEMIX + +C include Langmuir circulation parameterization +#undef ALLOW_GGL90_LANGMUIR + +#endif /* ALLOW_GGL90 */ +#endif /* GGL90_OPTIONS_H */ diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/GMREDI_OPTIONS.h b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/GMREDI_OPTIONS.h new file mode 100755 index 0000000..4ff989d --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/GMREDI_OPTIONS.h @@ -0,0 +1,69 @@ +#ifndef GMREDI_OPTIONS_H +#define GMREDI_OPTIONS_H +#include "PACKAGES_CONFIG.h" +#include "CPP_OPTIONS.h" + +CBOP +C !ROUTINE: GMREDI_OPTIONS.h +C !INTERFACE: +C #include "GMREDI_OPTIONS.h" + +C !DESCRIPTION: +C *==================================================================* +C | CPP options file for GM/Redi package: +C | Control which optional features to compile in this package code. +C *==================================================================* +CEOP + +#ifdef ALLOW_GMREDI +C Package-specific Options & Macros go here + +C Designed to simplify the Ajoint code: +#define GMREDI_WITH_STABLE_ADJOINT +C -- exclude the clipping/tapering part of the code that is not used +#define GM_EXCLUDE_CLIPPING +#define GM_EXCLUDE_FM07_TAP +#define GM_EXCLUDE_AC02_TAP +C #define GM_EXCLUDE_TAPERING +#define GM_EXCLUDE_SUBMESO + +C Allows to read-in background 3-D Redi and GM diffusivity coefficients +C Note: need these to be defined for use as control (pkg/ctrl) parameters +#define GM_READ_K3D_REDI +#define GM_READ_K3D_GM + +C This allows to use Visbeck et al formulation to compute K_GM+Redi +#undef GM_VISBECK_VARIABLE_K +C Use old calculation (before 2007/05/24) of Visbeck etal K_GM+Redi +C (which depends on tapering scheme) +#undef OLD_VISBECK_CALC + +C This allows the Bates et al formulation to calculate the +C bolus transport and K for Redi +#undef GM_BATES_K3D +#undef GM_BATES_PASSIVE + +C This allows the leading diagonal (top two rows) to be non-unity +C (a feature required when tapering adiabatically). +#define GM_NON_UNITY_DIAGONAL + +C Allows to use different values of K_GM and K_Redi ; also to +C be used with the advective form (Bolus velocity) of GM +#define GM_EXTRA_DIAGONAL + +C Allows to use the advective form (Bolus velocity) of GM +C instead of the Skew-Flux form (=default) +#define GM_BOLUS_ADVEC + +C Allows to use the Boundary-Value-Problem method to evaluate GM Bolus transport +#undef GM_BOLUS_BVP + +C Allow QG Leith variable viscosity to be added to GMRedi coefficient +#undef ALLOW_GM_LEITH_QG + +C Related to Adjoint-code: +#undef GM_AUTODIFF_EXCESSIVE_STORE +#undef GMREDI_MASK_SLOPES + +#endif /* ALLOW_GMREDI */ +#endif /* GMREDI_OPTIONS_H */ diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/MOM_COMMON_OPTIONS.h b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/MOM_COMMON_OPTIONS.h new file mode 100755 index 0000000..f80a356 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/MOM_COMMON_OPTIONS.h @@ -0,0 +1,25 @@ +C CPP options file for mom_common package +C Use this file for selecting CPP options within the mom_common package + +#ifndef MOM_COMMON_OPTIONS_H +#define MOM_COMMON_OPTIONS_H +#include "PACKAGES_CONFIG.h" +#include "CPP_OPTIONS.h" + +#ifdef ALLOW_MOM_COMMON +C Package-specific options go here + +C allow LeithQG coefficient to be calculated +#undef ALLOW_LEITH_QG + +C allow isotropic 3-D Smagorinsky viscosity +#undef ALLOW_SMAG_3D + +C allow full 3D specification of horizontal Laplacian Viscosity +#define ALLOW_3D_VISCAH + +C allow full 3D specification of horizontal Biharmonic Viscosity +#define ALLOW_3D_VISCA4 + +#endif /* ALLOW_MOM_COMMON */ +#endif /* MOM_COMMON_OPTIONS_H */ diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/SEAICE_LAYERS.h b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/SEAICE_LAYERS.h new file mode 100755 index 0000000..648c7b1 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/SEAICE_LAYERS.h @@ -0,0 +1,86 @@ +CBOP +C !ROUTINE: SEAICE_LAYERS.h + +C !DESCRIPTION: \bv +C *==========================================================* +C | SEAICE_LAYERS.h +C | o header file for sea ice multi-layer variables +C | for refined thermodynamics with vertical discretisation +C | of seaice and snow cover, accounting for heat content +C | and brine pockets contribution to seaice enthalpy +C *==========================================================* +C | Note: for now, only used when coupled to to GEOS AGCM +C *==========================================================* +C \ev +CEOP + +#ifdef HACK_FOR_GMAO_CPL +C SIqIce :: Seaice enthalpy for each ice layer and each category [J/m^2] +C SIqSnow :: Snow enthalpy for each snow layer and each category [J/m^2] +C SImeltPd :: Melt Pond volume for each category [m] +C SIiceAge :: Seaice Age for each category [s] +C SIskinS :: seaice skin salinity [psu] +C SIskinH :: seaice skin-layer depth [m] + _RL SIqIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nIceLayers, + & nITD,nSx,nSy) + _RL SIqSnow (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSnowLayers, + & nITD,nSx,nSy) + _RL SImeltPd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nITD,nSx,nSy) + _RL SIiceAge(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nITD,nSx,nSy) + _RL SIskinS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL SIskinH (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + COMMON /SEAICE_LAYERS/ + & SIqIce, SIqSnow, + & SImeltPd, SIiceAge, + & SIskinS, SIskinH + +C SIwindTauX :: wind stress over seaice, X-component at U point +C SIwindTauY :: wind stress over seaice, Y-component at V point + _RL SIwindTauX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL SIwindTauY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + COMMON /SEAICE_DYN_FORCING/ + & SIwindTauX, SIwindTauY + +C SIadv_Area :: advection increment of Seaice fraction [-] +C SIadv_Heff :: advection increment of Seaice thickness [m] +C SIadv_Hsnow :: advection increment of snow thickness [m] +C SIadv_tIces :: advection increment of ice surface temperature +C SIadv_qIce :: advection increment of Seaice enthalpy [J/m^2] +C SIadv_qSnow :: advection increment of Snow enthalpy [J/m^2] +C SIadv_meltPd :: advection increment of Melt Pond volume [m] +C SIadv_iceAge :: advection increment of Seaice Age [s] +C SIadv_skinS :: advection increment of seaice skin salinity [psu] +C SIadv_skinH :: advection increment of seaice skin-layer depth [m] + COMMON /SEAICE_ADV_INCREMENT/ + & SIadv_Area, SIadv_Heff, SIadv_Hsnow, + & SIadv_tIces, SIadv_qIce, SIadv_qSnow, + & SIadv_meltPd, SIadv_iceAge, SIadv_skinS, SIadv_skinH + _RL SIadv_Area (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nITD,nSx,nSy) + _RL SIadv_Heff (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nITD,nSx,nSy) + _RL SIadv_Hsnow (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nITD,nSx,nSy) + _RL SIadv_tIces (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nITD,nSx,nSy) + _RL SIadv_qIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nIceLayers, + & nITD,nSx,nSy) + _RL SIadv_qSnow (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSnowLayers, + & nITD,nSx,nSy) + _RL SIadv_meltPd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nITD,nSx,nSy) + _RL SIadv_iceAge(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nITD,nSx,nSy) + _RL SIadv_skinS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL SIadv_skinH (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + +C oceWeight :: grid-cell ocean fraction from GEOS [0-1] + _RL oceWeight(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + COMMON /DIAGS_GMAO_CPL/ + & oceWeight + +C SI_FRZMLT :: available heat (W/m^2) to freeze (>0) or melt (<0) sea ice +C so that surface level ocean reaches freezing temperature + _RL SI_FRZMLT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + COMMON /SEAICE_FRZMLT/ + & SI_FRZMLT + +#endif /* HACK_FOR_GMAO_CPL */ + +CEH3 ;;; Local Variables: *** +CEH3 ;;; mode:fortran *** +CEH3 ;;; End: *** diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/SEAICE_OPTIONS.h b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/SEAICE_OPTIONS.h new file mode 100755 index 0000000..71ae446 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/SEAICE_OPTIONS.h @@ -0,0 +1,264 @@ +#ifndef SEAICE_OPTIONS_H +#define SEAICE_OPTIONS_H +#include "PACKAGES_CONFIG.h" +#include "CPP_OPTIONS.h" + +C *==========================================================* +C | SEAICE_OPTIONS.h +C | o CPP options file for sea ice package. +C *==========================================================* +C | Use this file for selecting options within the sea ice +C | package. +C *==========================================================* + +#ifdef ALLOW_SEAICE +C--- Package-specific Options & Macros go here + +C-- Write "text-plots" of certain fields in STDOUT for debugging. +#undef SEAICE_DEBUG + +C-- By default, the sea-ice package uses its own integrated bulk +C formulae to compute fluxes (fu, fv, EmPmR, Qnet, and Qsw) over +C open-ocean. When this flag is set, these variables are computed +C in a separate external package, for example, pkg/exf, and then +C modified for sea-ice effects by pkg/seaice. +#define SEAICE_EXTERNAL_FLUXES + +C-- This CPP flag has been retired. The number of ice categories +C used to solve for seaice flux is now specified by run-time +C parameter SEAICE_multDim. +C Note: be aware of pickup_seaice.* compatibility issues when +C restarting a simulation with a different number of categories. +c#define SEAICE_MULTICATEGORY + +C-- run with sea Ice Thickness Distribution (ITD); +C set number of categories (nITD) in SEAICE_SIZE.h +#define SEAICE_ITD + +C-- Since the missing sublimation term is now included +C this flag is needed for backward compatibility +#undef SEAICE_DISABLE_SUBLIM + +C-- Suspected missing term in coupled ocn-ice heat budget (to be confirmed) +#undef SEAICE_DISABLE_HEATCONSFIX + +C-- Default is constant seaice salinity (SEAICE_salt0); Define the following +C flag to consider (space & time) variable salinity: advected and forming +C seaice with a fraction (=SEAICE_saltFrac) of freezing seawater salinity. +C- Note: SItracer also offers an alternative way to handle variable salinity. +#undef SEAICE_VARIABLE_SALINITY + +C-- Enable grease ice parameterization (requires to define ALLOW_SITRACER): +C The grease ice parameterization delays formation of solid sea ice from +C frazil ice by a time constant and provides a dynamic calculation of the +C initial solid sea ice thickness HO as a function of winds, currents and +C available grease ice volume. Grease ice does not significantly reduce heat +C loss from the ocean in winter and area covered by grease is thus handled +C like open water (For details see Smedsrud and Martin, 2014, Ann.Glac.). +C Set SItrName(1) = 'grease' in namelist SEAICE_PARM03 in data.seaice +C then output SItr01 is SItrNameLong(1) = 'grease ice volume fraction', +C with SItrUnit(1) = '[0-1]', which needs to be multiplied by SIheff to +C yield grease ice volume. Additionally, the actual grease ice layer +C thickness (diagnostic SIgrsLT) can be saved. +#undef SEAICE_GREASE + +C-- Tracers of ice and/or ice cover. +#ifdef SEAICE_GREASE +C SEAICE_GREASE code requires to define ALLOW_SITRACER +# define ALLOW_SITRACER +#else +# undef ALLOW_SITRACER +#endif +#ifdef ALLOW_SITRACER +C- To try avoid 'spontaneous generation' of tracer maxima by advdiff. +# define ALLOW_SITRACER_ADVCAP + +C- Include code to diagnose sea ice tracer budgets in +C seaice_advdiff.F and seaice_tracer_phys.F. Diagnostics are +C computed the "call diagnostics_fill" statement is commented out. +# undef ALLOW_SITRACER_DEBUG_DIAG +#endif /* ALLOW_SITRACER */ + +C-- Historically, the seaice model was discretized on a B-Grid. This +C discretization should still work but it is not longer actively tested +C and supported. The following flag should always be set in order to use +C the operational C-grid discretization. +#define SEAICE_CGRID + +#ifdef SEAICE_CGRID +C-- Options for the C-grid version only: + +C enable advection of sea ice momentum +# undef SEAICE_ALLOW_MOM_ADVECTION + +C enable JFNK code by defining the following flag +# undef SEAICE_ALLOW_JFNK +C enable Krylov code by defining the following flag +# undef SEAICE_ALLOW_KRYLOV + +C-- Use a different order when mapping 2D velocity arrays to 1D vector +C before passing it to FGMRES. +# undef SEAICE_JFNK_MAP_REORDER + +C to reproduce old verification results for JFNK +# undef SEAICE_PRECOND_EXTRA_EXCHANGE + +C enable LSR to use global (multi-tile) tri-diagonal solver +# undef SEAICE_GLOBAL_3DIAG_SOLVER + +C enable EVP code by defining the following flag +# define SEAICE_ALLOW_EVP +# ifdef SEAICE_ALLOW_EVP +C- When set use SEAICE_zetaMin and SEAICE_evpDampC to limit viscosities +C from below and above in seaice_evp: not necessary, and not recommended +# undef SEAICE_ALLOW_CLIPZETA + +C Include code to avoid underflows in EVP-code (copied from CICE). +C Many compilers can handle this more efficiently with the help of a flag. +# undef SEAICE_EVP_ELIMINATE_UNDERFLOWS + +C Include code to print residual of EVP iteration for debugging/diagnostics +# undef ALLOW_SEAICE_EVP_RESIDUAL +# endif /* SEAICE_ALLOW_EVP */ + +C smooth regularization (without max-function) of delta for +C better differentiability +# undef SEAICE_DELTA_SMOOTHREG + +C regularize zeta to zmax with a smooth tanh-function instead +C of a min(zeta,zmax). This improves convergence of iterative +C solvers (Lemieux and Tremblay 2009, JGR). No effect on EVP +# undef SEAICE_ZETA_SMOOTHREG + +C-- Different yield curves within the VP rheology framework +C allow the truncated ellipse rheology (runtime flag SEAICEuseTEM) +# undef SEAICE_ALLOW_TEM + +C allow the use of the Mohr Coulomb rheology (runtime flag SEAICEuseMCS) +C as defined in (Ip 1991) /!\ This is known to give unstable results, +C use with caution +# undef SEAICE_ALLOW_MCS + +C allow the use of Mohr Coulomb with elliptical plastic potential +C (runtime flag SEAICEuseMCE) +# undef SEAICE_ALLOW_MCE + +C allow the teardrop and parabolic lens rheology +C (runtime flag SEAICEuseTD and SEAICEusePL) +# undef SEAICE_ALLOW_TEARDROP + +C-- LSR solver settings +C Use LSR vector code; not useful on non-vector machines, because it +C slows down convergence considerably, but the extra iterations are +C more than made up by the much faster code on vector machines. For +C the only regularly test vector machine these flags a specified +C in the build options file SUPER-UX_SX-8_sxf90_awi, so that we comment +C them out here. +# undef SEAICE_VECTORIZE_LSR + +C Use zebra-method (alternate lines) for line-successive-relaxation +C This modification improves the convergence of the vector code +C dramatically, so that is may actually be useful in general, but +C that needs to be tested. Can be used without vectorization options. +# undef SEAICE_LSR_ZEBRA + +C Include code to print residual of nonlinear outer loop of LSR +# undef SEAICE_ALLOW_CHECK_LSR_CONVERGENCE + +C This flag is also required for an actual adjoint of seaice_lsr; +C increases memory requirements a lot. +# undef SEAICE_LSR_ADJOINT_ITER + +C Use parameterisation of grounding ice for a better representation +C of fastice in shallow seas +# undef SEAICE_ALLOW_BOTTOMDRAG + +#else /* not SEAICE_CGRID, but old B-grid */ +C-- Options for the B-grid version only: + +C- By default for B-grid dynamics solver wind stress under sea-ice is +C set to the same value as it would be if there was no sea-ice. +C Define following CPP flag for B-grid ice-ocean stress coupling. +# define SEAICE_BICE_STRESS + +C- By default for B-grid dynamics solver surface tilt is obtained +C indirectly via geostrophic velocities. Define following CPP +C in order to use ETAN instead. +# define EXPLICIT_SSH_SLOPE + +C- Defining this flag turns on FV-discretization of the B-grid LSOR solver. +C It is smoother and includes all metric terms, similar to C-grid solvers. +C It is here for completeness, but its usefulness is unclear. +# undef SEAICE_LSRBNEW + +#endif /* SEAICE_CGRID */ + +C-- Some regularisations +C- When set limit the Ice-Loading to mass of 1/5 of Surface ocean grid-box +#undef SEAICE_CAP_ICELOAD + +C- When set use SEAICE_clipVelocties = .true., to clip U/VICE at 40cm/s, +C not recommended +#undef SEAICE_ALLOW_CLIPVELS + +C- When set cap the sublimation latent heat flux in solve4temp according +C to the available amount of ice+snow. Otherwise this term is treated +C like all of the others -- residuals heat and fw stocks are passed to +C the ocean at the end of seaice_growth in a conservative manner. +C SEAICE_CAP_SUBLIM is not needed as of now, but kept just in case. +#undef SEAICE_CAP_SUBLIM + +C-- AD flags +C- TAF related flag, currently only used in seaice_ad_check_lev[1-4]_dir.h; +C it is unclear if this is ever needed. +#undef AUTODIFF_SOMETIMES_NEEDED + +C- Reset fields to zero to stabilise AD code of dynamics solver +C (resulting in wrong gradients) +#undef SEAICE_DYN_STABLE_ADJOINT + +C- Another flag to simplify dependencies for TAF-generated AD-code +C the thermodynamic part, mostly by resetting variables to zero +#undef SEAICE_MODIFY_GROWTH_ADJ + +C- Special seaice flag for AD testing +#undef SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING + +C-- Use the adjointable sea-ice thermodynamic model +C in seaice_growth_adx.F instead of seaice_growth.F +#undef SEAICE_USE_GROWTH_ADX + +C-- These flags are not strictly AD-related but may help obtaining +C simpler AD-code: +C- Do not compile code that resets AREA (or AREAITD) to a mininum value +C of SEAICE_area_floor (=SIeps with default of 1e-5) if there is +C some finite sea ice thickness +#undef DISABLE_AREA_FLOOR + +C- Do not compile growth/thermodynamics code (avoiding this code can +C also be done by setting runtime parameter usePWthermodynamics=F) +#undef DISABLE_SEAICE_GROWTH + +C- Allow sea-ice dynamic code. This option is provided so that, +C if turned off (#undef), to compile (and process with TAF) only the +C the thermodynamics component of the code. Note that, if needed, +C sea-ice dynamics can be turned off at runtime (SEAICEuseDYNAMICS=F). +#define SEAICE_ALLOW_DYNAMICS + +C- Do not compile/use seaice-related obcs code when using obcs. +#undef DISABLE_SEAICE_OBCS + +C-- Enable free drift code +#undef SEAICE_ALLOW_FREEDRIFT + +C-- pkg/seaice cost functions compile flags +C- Sea-ice volume (requires pkg/cost) +#undef ALLOW_COST_ICE +#ifdef ALLOW_COST_ICE +C- Enable template for sea-ice volume export in seaice_cost_export.F +C (requires pkg/cost & ALLOW_COST_ICE defined) +# undef ALLOW_SEAICE_COST_EXPORT +#endif /* ALLOW_COST_ICE */ + +#endif /* ALLOW_SEAICE */ +#endif /* SEAICE_OPTIONS_H */ diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/SEAICE_SIZE.h b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/SEAICE_SIZE.h new file mode 100755 index 0000000..671ec1f --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/SEAICE_SIZE.h @@ -0,0 +1,54 @@ +#ifdef ALLOW_SEAICE + +CBOP +C !ROUTINE: SEAICE_SIZE.h +C !INTERFACE: +C #include SEAICE_SIZE.h + +C !DESCRIPTION: +C Contains seaice array-size definition (number of tracers,categories). + +C SItrMaxNum :: number of passive tracers to allocate +C nITD :: number of seaice categories to allocate +CEOP + +C- Maximum Number of categories + INTEGER nITD +C-- +#ifdef SEAICE_ITD +C nITD defines number of ice thickness categories, +C i.e. size of additional dimension to AREA, HEFF, HSNOW, etc. +C Bitz et al. (2001, JGR) suggest a minimum of nITD = 5 + PARAMETER (nITD = 5) +#else + PARAMETER (nITD = 7) +#endif + +C- Maximum Number of tracers + INTEGER SItrMaxNum + PARAMETER(SItrMaxNum = 3 ) + +#ifdef HACK_FOR_GMAO_CPL +C nIceLayers :: number of Ice layers (in each category) +C nSnowLayers :: number of snow layers (in each category) + INTEGER nIceLayers, nSnowLayers + PARAMETER( nIceLayers = 4 , nSnowLayers = 1 ) +#endif /* HACK_FOR_GMAO_CPL */ + +#ifdef ALLOW_AUTODIFF + INTEGER iicekey + INTEGER nEVPstepMax + PARAMETER ( nEVPstepMax=180 ) + INTEGER NMAX_TICE + PARAMETER ( NMAX_TICE=10 ) + INTEGER SOLV_MAX_FIXED + PARAMETER ( SOLV_MAX_FIXED=500 ) + INTEGER MPSEUDOTIMESTEPS + PARAMETER (MPSEUDOTIMESTEPS=2) +#endif /* ALLOW_AUTODIFF */ + +#endif /* ALLOW_SEAICE */ + +CEH3 ;;; Local Variables: *** +CEH3 ;;; mode:fortran *** +CEH3 ;;; End: *** diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/SIZE.h b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/SIZE.h new file mode 100755 index 0000000..8f8593d --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/SIZE.h @@ -0,0 +1,65 @@ +CBOP +C !ROUTINE: SIZE.h +C !INTERFACE: +C include SIZE.h +C !DESCRIPTION: \bv +C *==========================================================* +C | SIZE.h Declare size of underlying computational grid. +C *==========================================================* +C | The design here supports a three-dimensional model grid +C | with indices I,J and K. The three-dimensional domain +C | is comprised of nPx*nSx blocks (or tiles) of size sNx +C | along the first (left-most index) axis, nPy*nSy blocks +C | of size sNy along the second axis and one block of size +C | Nr along the vertical (third) axis. +C | Blocks/tiles have overlap regions of size OLx and OLy +C | along the dimensions that are subdivided. +C *==========================================================* +C \ev +C +C Voodoo numbers controlling data layout: +C sNx :: Number of X points in tile. +C sNy :: Number of Y points in tile. +C OLx :: Tile overlap extent in X. +C OLy :: Tile overlap extent in Y. +C nSx :: Number of tiles per process in X. +C nSy :: Number of tiles per process in Y. +C nPx :: Number of processes to use in X. +C nPy :: Number of processes to use in Y. +C Nx :: Number of points in X for the full domain. +C Ny :: Number of points in Y for the full domain. +C Nr :: Number of points in vertical direction. +CEOP + INTEGER sNx + INTEGER sNy + INTEGER OLx + INTEGER OLy + INTEGER nSx + INTEGER nSy + INTEGER nPx + INTEGER nPy + INTEGER Nx + INTEGER Ny + INTEGER Nr + PARAMETER ( + & sNx = 30, + & sNy = 30, + & OLx = 4, + & OLy = 4, + & nSx = 1, + & nSy = 1, +! & nPx = 1053-286, + & nPx = 768, + & nPy = 1, + & Nx = sNx*nSx*nPx, + & Ny = sNy*nSy*nPy, + & Nr = 50 ) + +C MAX_OLX :: Set to the maximum overlap region size of any array +C MAX_OLY that will be exchanged. Controls the sizing of exch +C routine buffers. + INTEGER MAX_OLX + INTEGER MAX_OLY + PARAMETER ( MAX_OLX = OLx, + & MAX_OLY = OLy ) + diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/mom_calc_visc.F b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/mom_calc_visc.F new file mode 100755 index 0000000..a9bda8e --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/mom_calc_visc.F @@ -0,0 +1,794 @@ +#include "MOM_COMMON_OPTIONS.h" +#ifdef ALLOW_AUTODIFF +# include "AUTODIFF_OPTIONS.h" +#endif + +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| +CBOP +C !ROUTINE: MOM_CALC_VISC + +C !INTERFACE: + SUBROUTINE MOM_CALC_VISC( + I bi,bj,k, + O viscAh_Z,viscAh_D,viscA4_Z,viscA4_D, + I hDiv,vort3,tension,strain,stretching,KE,hFacZ, + I myThid) + +C !DESCRIPTION: +C Calculate horizontal viscosities (L is typical grid width) +C harmonic viscosity= +C viscAh (or viscAhD on div pts and viscAhZ on zeta pts) +C +0.25*L**2*viscAhGrid/deltaT +C +sqrt((viscC2leith/pi)**6*grad(Vort3)**2 +C +(viscC2leithD/pi)**6*grad(hDiv)**2)*L**3 +C +(viscC2smag/pi)**2*L**2*sqrt(Tension**2+Strain**2) +C +C biharmonic viscosity= +C viscA4 (or viscA4D on div pts and viscA4Z on zeta pts) +C +0.25*0.125*L**4*viscA4Grid/deltaT (approx) +C +0.125*L**5*sqrt((viscC4leith/pi)**6*grad(Vort3)**2 +C +(viscC4leithD/pi)**6*grad(hDiv)**2) +C +0.125*L**4*(viscC4smag/pi)**2*sqrt(Tension**2+Strain**2) +C +C Note that often 0.125*L**2 is the scale between harmonic and +C biharmonic (see Griffies and Hallberg (2000)) +C This allows the same value of the coefficient to be used +C for roughly similar results with biharmonic and harmonic +C +C LIMITERS -- limit min and max values of viscosities +C viscAhReMax is min value for grid point harmonic Reynolds num +C harmonic viscosity>sqrt(2*KE)*L/viscAhReMax +C +C viscA4ReMax is min value for grid point biharmonic Reynolds num +C biharmonic viscosity>sqrt(2*KE)*L**3/8/viscA4ReMax +C +C viscAhgridmax is CFL stability limiter for harmonic viscosity +C harmonic viscosity<0.25*viscAhgridmax*L**2/deltaT +C +C viscA4gridmax is CFL stability limiter for biharmonic viscosity +C biharmonic viscosity0.25*viscAhgridmin*L**2/deltaT +C biharmonic viscosity>viscA4gridmin*L**4/32/deltaT (approx) + +C RECOMMENDED VALUES +C viscC2Leith=1-3 +C viscC2LeithD=1-3 +C viscC2LeithQG=1 +C viscC4Leith=1-3 +C viscC4LeithD=1.5-3 +C viscC2smag=2.2-4 (Griffies and Hallberg,2000) +C 0.2-0.9 (Smagorinsky,1993) +C viscC4smag=2.2-4 (Griffies and Hallberg,2000) +C viscAhReMax>=1, (<2 suppresses a computational mode) +C viscA4ReMax>=1, (<2 suppresses a computational mode) +C viscAhgridmax=1 +C viscA4gridmax=1 +C viscAhgrid<1 +C viscA4grid<1 +C viscAhgridmin<<1 +C viscA4gridmin<<1 + +C !USES: + IMPLICIT NONE + +C == Global variables == +#include "SIZE.h" +#include "GRID.h" +#include "EEPARAMS.h" +#include "PARAMS.h" +#include "MOM_VISC.h" +#ifdef ALLOW_AUTODIFF_TAMC +#include "tamc.h" +#endif /* ALLOW_AUTODIFF_TAMC */ + +C !INPUT/OUTPUT PARAMETERS: +C myThid :: my thread Id number + INTEGER bi,bj,k + _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscA4_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscA4_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL tension(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL strain(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL stretching(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + INTEGER myThid +CEOP + +C !LOCAL VARIABLES: + INTEGER i,j +#ifdef ALLOW_NONHYDROSTATIC + _RL shiftAh, shiftA4 +#endif +#ifdef ALLOW_AUTODIFF_TAMC +C kkey :: tape key (depends on levels and tiles) +C ijkkey :: tape key (depends on i,j-indices, levels, and tiles) + INTEGER kkey, ijkkey +#endif + _RL smag2fac, smag4fac + _RL leith2fac, leith4fac + _RL leithD2fac, leithD4fac + _RL leithQG2fac + _RL viscAhRe_max, viscA4Re_max + _RL Alin, grdVrt, grdDiv, keZpt + _RL deepFac3, deepFac4 + _RL L2, L3, L5, L2rdt, L4rdt, recip_dt + _RL Uscl,U4scl + _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL vrtDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL vrtDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscA4_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscAh_ZMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscAh_DMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscA4_ZMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscA4_DMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscAh_ZLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscAh_DLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscA4_ZLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscA4_DLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscAh_ZLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscAh_DLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscA4_ZLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscA4_DLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy) +#ifdef ALLOW_LEITH_QG + _RL viscAh_ZLthQG(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscAh_DLthQG(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL sqargQG +#endif + _RL viscAh_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscAh_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscA4_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL viscA4_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL sqargAh, sqargA4, sqargAhD, sqargA4D, sqargSmag + LOGICAL calcLeith, calcSmag, calcLeithQG + +#ifdef ALLOW_AUTODIFF_TAMC + kkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy + kkey = k + (kkey-1)*Nr +#endif /* ALLOW_AUTODIFF_TAMC */ + +C-- Set flags which are used in this S/R and elsewhere : +C useVariableVisc, useHarmonicVisc and useBiharmonicVisc +C are now set early on (in S/R SET_PARAMS) + +c IF ( useVariableVisc ) THEN +C---- variable viscosity : + + recip_dt = 1. _d 0 + IF ( deltaTMom.NE.0. ) recip_dt = 1. _d 0/deltaTMom + deepFac3 = deepFac2C(k)*deepFacC(k) + deepFac4 = deepFac2C(k)*deepFac2C(k) + + IF ( useHarmonicVisc .AND. viscAhReMax.NE.0. ) THEN + viscAhRe_max=SQRT(2. _d 0)/viscAhReMax + ELSE + viscAhRe_max=0. _d 0 + ENDIF + + IF ( useBiharmonicVisc .AND. viscA4ReMax.NE.0. ) THEN + viscA4Re_max=0.125 _d 0*SQRT(2. _d 0)/viscA4ReMax + ELSE + viscA4Re_max=0. _d 0 + ENDIF + + calcLeithQG = (viscC2LeithQG.NE.zeroRL) + calcLeith= + & (viscC2leith.NE.0.) + & .OR.(viscC2leithD.NE.0.) + & .OR.(viscC4leith.NE.0.) + & .OR.(viscC4leithD.NE.0.) + & .OR. calcLeithQG + + calcSmag= + & (viscC2smag.NE.0.) + & .OR.(viscC4smag.NE.0.) + + IF (calcSmag) THEN + smag2fac=(viscC2smag/pi)**2 + smag4fac=0.125 _d 0*(viscC4smag/pi)**2 + ELSE + smag2fac=0. _d 0 + smag4fac=0. _d 0 + ENDIF + + IF (calcLeith) THEN + IF (useFullLeith) THEN +C Uses correct calculation for gradients, but might not work on cube sphere + leith2fac =(viscC2leith /pi)**6 + leithD2fac=(viscC2leithD/pi)**6 + leithQG2fac = (viscC2LeithQG/pi)**6 + leith4fac =0.015625 _d 0*(viscC4leith /pi)**6 + leithD4fac=0.015625 _d 0*(viscC4leithD/pi)**6 + ELSE +C Uses approximate gradients, but works on cube sphere. No reason to use this +C unless `useFullLeith` fails for your setup. + leith2fac =(viscC2leith /pi)**3 + leithD2fac=(viscC2leithD/pi)**3 + leithQG2fac = (viscC2LeithQG/pi)**3 + leith4fac =0.125 _d 0*(viscC4leith /pi)**3 + leithD4fac=0.125 _d 0*(viscC4leithD/pi)**3 + ENDIF + ELSE + leith2fac=0. _d 0 + leith4fac=0. _d 0 + leithQG2fac=0. _d 0 + leithD2fac=0. _d 0 + leithD4fac=0. _d 0 + ENDIF + + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx +C- viscosity arrays have been initialised everywhere before calling this S/R +c viscAh_D(i,j) = viscAhD +c viscAh_Z(i,j) = viscAhZ +c viscA4_D(i,j) = viscA4D +c viscA4_Z(i,j) = viscA4Z + + viscAh_DLth(i,j) = 0. _d 0 + viscAh_ZLth(i,j) = 0. _d 0 + viscA4_DLth(i,j) = 0. _d 0 + viscA4_ZLth(i,j) = 0. _d 0 + viscAh_DLthD(i,j)= 0. _d 0 + viscAh_ZLthD(i,j)= 0. _d 0 + viscA4_DLthD(i,j)= 0. _d 0 + viscA4_ZLthD(i,j)= 0. _d 0 +#ifdef ALLOW_LEITH_QG + viscAh_DLthQG(i,j) = 0. _d 0 + viscAh_ZLthQG(i,j) = 0. _d 0 +#endif + + viscAh_DSmg(i,j) = 0. _d 0 + viscAh_ZSmg(i,j) = 0. _d 0 + viscA4_DSmg(i,j) = 0. _d 0 + viscA4_ZSmg(i,j) = 0. _d 0 + ENDDO + ENDDO + +C- Initialise to zero gradient of vorticity and divergence: + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + divDx(i,j) = 0. + divDy(i,j) = 0. + vrtDx(i,j) = 0. + vrtDy(i,j) = 0. + ENDDO + ENDDO + + IF ( calcLeith ) THEN +C-- horizontal gradient of horizontal divergence: +C- gradient in x direction: + IF (useCubedSphereExchange) THEN +C to compute d/dx(hDiv), fill corners with appropriate values: + CALL FILL_CS_CORNER_TR_RL( 1, .FALSE., + & hDiv, bi,bj, myThid ) + ENDIF + DO j=2-OLy,sNy+OLy-1 + DO i=2-OLx,sNx+OLx-1 + divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j)) + & *recip_dxC(i,j,bi,bj)*recip_deepFacC(k) + ENDDO + ENDDO + +C- gradient in y direction: + IF (useCubedSphereExchange) THEN +C to compute d/dy(hDiv), fill corners with appropriate values: + CALL FILL_CS_CORNER_TR_RL( 2, .FALSE., + & hDiv, bi,bj, myThid ) + ENDIF + DO j=2-OLy,sNy+OLy-1 + DO i=2-OLx,sNx+OLx-1 + divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1)) + & *recip_dyC(i,j,bi,bj)*recip_deepFacC(k) + ENDDO + ENDDO + +C-- horizontal gradient of vertical vorticity: +C- gradient in x direction: + DO j=2-OLy,sNy+OLy + DO i=2-OLx,sNx+OLx-1 + vrtDx(i,j) = (vort3(i+1,j)-vort3(i,j)) + & *recip_dxG(i,j,bi,bj)*recip_deepFacC(k) + & *maskS(i,j,k,bi,bj) +#ifdef ALLOW_OBCS + & *maskInS(i,j,bi,bj) +#endif + ENDDO + ENDDO +C- gradient in y direction: + DO j=2-OLy,sNy+OLy-1 + DO i=2-OLx,sNx+OLx + vrtDy(i,j) = (vort3(i,j+1)-vort3(i,j)) + & *recip_dyG(i,j,bi,bj)*recip_deepFacC(k) + & *maskW(i,j,k,bi,bj) +#ifdef ALLOW_OBCS + & *maskInW(i,j,bi,bj) +#endif + ENDDO + ENDDO + +#ifdef ALLOW_LEITH_QG + IF ( calcLeithQG ) THEN +C horizontal gradient of vorticity and vortex stretching: +C In the case of using QG Leith, we want to add a term +C before calculating vector magnitude, so add to the +C values just calculated. +C gradient in x direction: + DO j=2-OLy,sNy+OLy + DO i=2-OLx,sNx+OLx-1 +C Average d/dx of stretching onto V-points to match vrtDX + vrtDx(i,j) = vrtDx(i,j) + & + halfRL*halfRL* + & ( (stretching(i+1,j)-stretching(i,j)) + & *recip_dxC(i+1,j,bi,bj)*recip_deepFacC(k) + & + (stretching(i,j)-stretching(i-1,j)) + & *recip_dxC(i,j,bi,bj)*recip_deepFacC(k) + & + (stretching(i+1,j-1)-stretching(i,j-1)) + & *recip_dxC(i,j-1,bi,bj)*recip_deepFacC(k) + & + (stretching(i,j-1)-stretching(i-1,j-1)) + & *recip_dxC(i-1,j-1,bi,bj)*recip_deepFacC(k) + & )*maskS(i,j,k,bi,bj) +#ifdef ALLOW_OBCS + & *maskInS(i,j,bi,bj) +#endif + ENDDO + ENDDO +C- gradient in y direction: + DO j=2-OLy,sNy+OLy-1 + DO i=2-OLx,sNx+OLx +C Average d/dy of stretching onto U-points to match vrtDy + vrtDy(i,j) = vrtDy(i,j) + & + halfRL*halfRL* + & ( (stretching(i,j+1)-stretching(i,j)) + & *recip_dyC(i,j+1,bi,bj)*recip_deepFacC(k) + & + (stretching(i,j)-stretching(i,j-1)) + & *recip_dyC(i,j,bi,bj)*recip_deepFacC(k) + & + (stretching(i-1,j+1)-stretching(i-1,j)) + & *recip_dyC(i-1,j+1,bi,bj)*recip_deepFacC(k) + & + (stretching(i-1,j)-stretching(i-1,j-1)) + & *recip_dyC(i-1,j,bi,bj)*recip_deepFacC(k) + & )*maskW(i,j,k,bi,bj) +#ifdef ALLOW_OBCS + & *maskInW(i,j,bi,bj) +#endif + ENDDO + ENDDO +C end if calcLeithQG + ENDIF +#endif /* ALLOW_LEITH_QG */ + +C-- end if calcLeith + ENDIF + + DO j=2-OLy,sNy+OLy-1 + DO i=2-OLx,sNx+OLx-1 +CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC + +#ifdef ALLOW_AUTODIFF_TAMC +# ifndef AUTODIFF_DISABLE_LEITH + ijkkey = i+OLx + (j+OLy-1)*(sNx+2*OLx) + & + (kkey-1)*(sNx+2*OLx)*(sNy+2*OLy) +CADJ STORE viscA4_ZSmg(i,j)=comlev1_mom_ijk_loop,key=ijkkey,byte=isbyte +CADJ STORE viscAh_ZSmg(i,j)=comlev1_mom_ijk_loop,key=ijkkey,byte=isbyte +# endif +#endif /* ALLOW_AUTODIFF_TAMC */ + +C These are (powers of) length scales + L2 = L2_D(i,j,bi,bj)*deepFac2C(k) + L2rdt = 0.25 _d 0*recip_dt*L2 + L3 = L3_D(i,j,bi,bj)*deepFac3 + L4rdt = L4rdt_D(i,j,bi,bj)*deepFac4 + L5 = (L2*L3) + +#ifndef AUTODIFF_DISABLE_REYNOLDS_SCALE +C Velocity Reynolds Scale + IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN + Uscl=SQRT(KE(i,j)*L2)*viscAhRe_max + ELSE + Uscl=0. + ENDIF + IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN + U4scl=SQRT(KE(i,j))*L3*viscA4Re_max + ELSE + U4scl=0. + ENDIF +#endif /* ndef AUTODIFF_DISABLE_REYNOLDS_SCALE */ + +#ifndef AUTODIFF_DISABLE_LEITH + IF (useFullLeith.AND.calcLeith) THEN +C This is the vector magnitude of the vorticity gradient squared + grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1) + & + vrtDx(i,j)*vrtDx(i,j) ) + & + (vrtDy(i+1,j)*vrtDy(i+1,j) + & + vrtDy(i,j)*vrtDy(i,j) ) ) + +C This is the vector magnitude of grad (div.v) squared +C Using it in Leith serves to damp instabilities in w. + grdDiv=0.25 _d 0*( (divDx(i+1,j)*divDx(i+1,j) + & + divDx(i,j)*divDx(i,j) ) + & + (divDy(i,j+1)*divDy(i,j+1) + & + divDy(i,j)*divDy(i,j) ) ) + + sqargAh = leith2fac*grdVrt+leithD2fac*grdDiv + sqargA4 = leith4fac*grdVrt+leithD4fac*grdDiv + sqargAhD = leithD2fac*grdDiv + sqargA4D = leithD4fac*grdDiv +#ifdef ALLOW_LEITH_QG + sqargQG = leithQG2fac*(grdVrt+grdDiv) +#endif + +#ifdef ALLOW_AUTODIFF +C Avoid derivative of SQRT(0) + IF (sqargAh .GT.0. _d 0) sqargAh = SQRT(sqargAh) + IF (sqargA4 .GT.0. _d 0) sqargA4 = SQRT(sqargA4) + IF (sqargAhD .GT.0. _d 0) sqargAhD = SQRT(sqargAhD) + IF (sqargA4D .GT.0. _d 0) sqargA4D = SQRT(sqargA4D) +# ifdef ALLOW_LEITH_QG + IF (sqargQG .GT.0. _d 0) sqargQG = SQRT(sqargQG) +# endif +#else /* ALLOW_AUTODIFF */ + sqargAh = SQRT(sqargAh) + sqargA4 = SQRT(sqargA4) + sqargAhD = SQRT(sqargAhD) + sqargA4D = SQRT(sqargA4D) +# ifdef ALLOW_LEITH_QG + sqargQG = SQRT(sqargQG) +# endif +#endif /* ALLOW_AUTODIFF */ + viscAh_DLth(i,j) = sqargAh * L3 + viscA4_DLth(i,j) = sqargA4 * L5 + viscAh_DLthd(i,j)= sqargAhD * L3 + viscA4_DLthd(i,j)= sqargA4D * L5 +#ifdef ALLOW_LEITH_QG + viscAh_DLthQG(i,j)=sqargQG * L3 +#endif + + ELSEIF (calcLeith) THEN +C but this approximation will work on cube (and differs by as much as 4X) + grdVrt=MAX( ABS(vrtDx(i,j+1)), ABS(vrtDx(i,j)) ) + grdVrt=MAX( grdVrt, ABS(vrtDy(i+1,j)) ) + grdVrt=MAX( grdVrt, ABS(vrtDy(i,j)) ) + +C This approximation is good to the same order as above... + grdDiv=MAX( ABS(divDx(i+1,j)), ABS(divDx(i,j)) ) + grdDiv=MAX( grdDiv, ABS(divDy(i,j+1)) ) + grdDiv=MAX( grdDiv, ABS(divDy(i,j)) ) + + viscAh_DLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3 + viscA4_DLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5 + viscAh_DLthD(i,j)=((leithD2fac*grdDiv))*L3 + viscA4_DLthD(i,j)=((leithD4fac*grdDiv))*L5 +#ifdef ALLOW_LEITH_QG + viscAh_DLthQG(i,j)=leithQG2fac*(grdVrt + grdDiv)*L3 +#endif + + ELSE + viscAh_DLth(i,j)=0. _d 0 + viscA4_DLth(i,j)=0. _d 0 + viscAh_DLthD(i,j)=0. _d 0 + viscA4_DLthD(i,j)=0. _d 0 +#ifdef ALLOW_LEITH_QG + viscAh_DLthQG(i,j)=0. _d 0 +#endif + ENDIF + + IF (calcSmag) THEN + sqargSmag = tension(i,j)**2 + & +0.25 _d 0*(strain(i+1, j )**2+strain( i ,j+1)**2 + & +strain(i , j )**2+strain(i+1,j+1)**2) +#ifdef ALLOW_AUTODIFF +C Avoid derivative of SQRT(0) + IF (sqargSmag.GT.0. _d 0) sqargSmag = SQRT(sqargSmag) +#else + sqargSmag = SQRT(sqargSmag) +#endif + viscAh_DSmg(i,j)=L2*sqargSmag + viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j) + viscAh_DSmg(i,j)=smag2fac*viscAh_DSmg(i,j) + ELSE + viscAh_DSmg(i,j)=0. _d 0 + viscA4_DSmg(i,j)=0. _d 0 + ENDIF +#endif /* AUTODIFF_DISABLE_LEITH */ + +C Harmonic on Div.u points + Alin=viscAhD+viscAhGrid*L2rdt + & + viscAh_DLth(i,j)+viscAh_DSmg(i,j) +#ifdef ALLOW_LEITH_QG + & + viscAh_DLthQG(i,j) +#endif +#ifdef ALLOW_3D_VISCAH + & + viscAhDfld(i,j,k,bi,bj) +# ifdef AUTODIFF_ALLOW_VISCFACADJ + & *viscFacAdj +# endif /* AUTODIFF_ALLOW_VISCFACADJ */ +#endif /* ALLOW_3D_VISCAH */ + viscAh_DMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl) + viscAh_D(i,j)=MAX(viscAh_DMin(i,j),Alin) + viscAh_DMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax) + viscAh_D(i,j)=MIN(viscAh_DMax(i,j),viscAh_D(i,j)) + + if ( (yC(i,j,bi,bj).GE.33.) .AND. + & (yC(i,j,bi,bj).LE.39.) .AND. + & (xC(i,j,bi,bj).GE.-7.) .AND. + & (xC(i,j,bi,bj).LE.-2.) + & ) then + viscAh_D(i,j)=10. _d 0 * viscAh_D(i,j) + endif + +C BiHarmonic on Div.u points + Alin=viscA4D+viscA4Grid*L4rdt + & + viscA4_DLth(i,j)+viscA4_DSmg(i,j) +#ifdef ALLOW_3D_VISCA4 + & + viscA4Dfld(i,j,k,bi,bj) +# ifdef AUTODIFF_ALLOW_VISCFACADJ + & *viscFacAdj +# endif /* AUTODIFF_ALLOW_VISCFACADJ */ +#endif /* ALLOW_3D_VISCA4 */ + viscA4_DMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl) + viscA4_D(i,j)=MAX(viscA4_DMin(i,j),Alin) + viscA4_DMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max) + viscA4_D(i,j)=MIN(viscA4_DMax(i,j),viscA4_D(i,j)) + +CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC +C These are (powers of) length scales + L2 = L2_Z(i,j,bi,bj)*deepFac2C(k) + L2rdt = 0.25 _d 0*recip_dt*L2 + L3 = L3_Z(i,j,bi,bj)*deepFac3 + L4rdt = L4rdt_Z(i,j,bi,bj)*deepFac4 + L5 = (L2*L3) + +#ifndef AUTODIFF_DISABLE_REYNOLDS_SCALE +C Velocity Reynolds Scale (Pb here at CS-grid corners !) + IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN + keZpt=0.25 _d 0*( (KE(i,j)+KE(i-1,j-1)) + & +(KE(i-1,j)+KE(i,j-1)) ) + IF ( keZpt.GT.0. ) THEN + Uscl = SQRT(keZpt*L2)*viscAhRe_max + U4scl= SQRT(keZpt)*L3*viscA4Re_max + ELSE + Uscl =0. + U4scl=0. + ENDIF + ELSE + Uscl =0. + U4scl=0. + ENDIF +#endif /* ndef AUTODIFF_DISABLE_REYNOLDS_SCALE */ + +#ifndef AUTODIFF_DISABLE_LEITH +C This is the vector magnitude of the vorticity gradient squared + IF (useFullLeith.AND.calcLeith) THEN + grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j) + & + vrtDx(i,j)*vrtDx(i,j) ) + & + (vrtDy(i,j-1)*vrtDy(i,j-1) + & + vrtDy(i,j)*vrtDy(i,j) ) ) + +C This is the vector magnitude of grad(div.v) squared + grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1) + & + divDx(i,j)*divDx(i,j) ) + & + (divDy(i-1,j)*divDy(i-1,j) + & + divDy(i,j)*divDy(i,j) ) ) + + sqargAh = leith2fac*grdVrt+leithD2fac*grdDiv + sqargA4 = leith4fac*grdVrt+leithD4fac*grdDiv + sqargAhD = leithD2fac*grdDiv + sqargA4D = leithD4fac*grdDiv +#ifdef ALLOW_LEITH_QG + sqargQG = leithQG2fac*(grdVrt+grdDiv) +#endif +#ifdef ALLOW_AUTODIFF +C Avoid derivative of SQRT(0) + IF (sqargAh .GT.0. _d 0) sqargAh = SQRT(sqargAh) + IF (sqargA4 .GT.0. _d 0) sqargA4 = SQRT(sqargA4) + IF (sqargAhD .GT.0. _d 0) sqargAhD = SQRT(sqargAhD) + IF (sqargA4D .GT.0. _d 0) sqargA4D = SQRT(sqargA4D) +# ifdef ALLOW_LEITH_QG + IF (sqargQG .GT.0. _d 0) sqargQG = SQRT(sqargQG) +# endif +#else /* ALLOW_AUTODIFF */ + sqargAh = SQRT(sqargAh) + sqargA4 = SQRT(sqargA4) + sqargAhD = SQRT(sqargAhD) + sqargA4D = SQRT(sqargA4D) +# ifdef ALLOW_LEITH_QG + sqargQG = SQRT(sqargQG) +# endif +#endif /* ALLOW_AUTODIFF */ + viscAh_ZLth(i,j) = sqargAh * L3 + viscA4_ZLth(i,j) = sqargA4 * L5 + viscAh_ZLthd(i,j)= sqargAhD * L3 + viscA4_ZLthd(i,j)= sqargA4D * L5 +#ifdef ALLOW_LEITH_QG + viscAh_ZLthQG(i,j)=sqargQG * L3 +#endif + + ELSEIF (calcLeith) THEN +C but this approximation will work on cube (and differs by 4X) + grdVrt=MAX( ABS(vrtDx(i-1,j)), ABS(vrtDx(i,j)) ) + grdVrt=MAX( grdVrt, ABS(vrtDy(i,j-1)) ) + grdVrt=MAX( grdVrt, ABS(vrtDy(i,j)) ) + + grdDiv=MAX( ABS(divDx(i,j)), ABS(divDx(i,j-1)) ) + grdDiv=MAX( grdDiv, ABS(divDy(i,j)) ) + grdDiv=MAX( grdDiv, ABS(divDy(i-1,j)) ) + + viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3 + viscA4_ZLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5 + viscAh_ZLthD(i,j)=(leithD2fac*grdDiv)*L3 + viscA4_ZLthD(i,j)=(leithD4fac*grdDiv)*L5 +#ifdef ALLOW_LEITH_QG + viscAh_ZLthQG(i,j)=leithQG2fac*(grdVrt + grdDiv)*L3 +#endif + ELSE + viscAh_ZLth(i,j)=0. _d 0 + viscA4_ZLth(i,j)=0. _d 0 + viscAh_ZLthD(i,j)=0. _d 0 + viscA4_ZLthD(i,j)=0. _d 0 +#ifdef ALLOW_LEITH_QG + viscAh_ZLthQG(i,j)=0. _d 0 +#endif + ENDIF + + IF (calcSmag) THEN + sqargSmag = strain(i,j)**2 + & +0.25 _d 0*(tension( i , j )**2+tension( i ,j-1)**2 + & +tension(i-1, j )**2+tension(i-1,j-1)**2) +#ifdef ALLOW_AUTODIFF +C Avoid derivative of SQRT(0) + IF (sqargSmag.GT.0. _d 0) sqargSmag = SQRT(sqargSmag) +#else + sqargSmag = SQRT(sqargSmag) +#endif + viscAh_ZSmg(i,j)=L2*sqargSmag + viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j) + viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j) + ENDIF +#endif /* AUTODIFF_DISABLE_LEITH */ + +C Harmonic on Zeta points + Alin=viscAhZ+viscAhGrid*L2rdt + & + viscAh_ZLth(i,j)+viscAh_ZSmg(i,j) +#ifdef ALLOW_LEITH_QG + & + viscAh_ZLthQG(i,j) +#endif +#ifdef ALLOW_3D_VISCAH + & + viscAhZfld(i,j,k,bi,bj) +# ifdef AUTODIFF_ALLOW_VISCFACADJ + & *viscFacAdj +# endif /* AUTODIFF_ALLOW_VISCFACADJ */ +#endif + viscAh_ZMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl) + viscAh_Z(i,j)=MAX(viscAh_ZMin(i,j),Alin) + viscAh_ZMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax) + viscAh_Z(i,j)=MIN(viscAh_ZMax(i,j),viscAh_Z(i,j)) + + if ( (yG(i,j,bi,bj).GE.33.) .AND. + & (yG(i,j,bi,bj).LE.39.) .AND. + & (xG(i,j,bi,bj).GE.-7.) .AND. + & (xG(i,j,bi,bj).LE.-2.) + & ) then + viscAh_Z(i,j)=10. _d 0 * viscAh_Z(i,j) + endif + +C BiHarmonic on Zeta points + Alin=viscA4Z+viscA4Grid*L4rdt + & + viscA4_ZLth(i,j)+viscA4_ZSmg(i,j) +#ifdef ALLOW_3D_VISCA4 + & + viscA4Zfld(i,j,k,bi,bj) +# ifdef AUTODIFF_ALLOW_VISCFACADJ + & *viscFacAdj +# endif /* AUTODIFF_ALLOW_VISCFACADJ */ +#endif + viscA4_ZMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl) + viscA4_Z(i,j)=MAX(viscA4_ZMin(i,j),Alin) + viscA4_ZMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max) + viscA4_Z(i,j)=MIN(viscA4_ZMax(i,j),viscA4_Z(i,j)) + ENDDO + ENDDO + +#ifdef ALLOW_NONHYDROSTATIC + IF ( nonHydrostatic ) THEN +C-- Pass Viscosities to calc_gw (if constant, not necessary) + + IF ( k.LT.Nr ) THEN +C Prepare for next level (next call) + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + viscAh_W(i,j,k+1,bi,bj) = halfRL*viscAh_D(i,j) + viscA4_W(i,j,k+1,bi,bj) = halfRL*viscA4_D(i,j) + ENDDO + ENDDO + ENDIF + + shiftAh = viscAhW - viscAhD + shiftA4 = viscA4W - viscA4D + IF ( k.EQ.1 ) THEN +C These values dont get used + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + viscAh_W(i,j,k,bi,bj) = shiftAh + viscAh_D(i,j) + viscA4_W(i,j,k,bi,bj) = shiftA4 + viscA4_D(i,j) + ENDDO + ENDDO + ELSE +C Note that previous call of this function has already added half. + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + viscAh_W(i,j,k,bi,bj) = shiftAh + viscAh_W(i,j,k,bi,bj) + & + halfRL*viscAh_D(i,j) + viscA4_W(i,j,k,bi,bj) = shiftA4 + viscA4_W(i,j,k,bi,bj) + & + halfRL*viscA4_D(i,j) + ENDDO + ENDDO + ENDIF + + ENDIF +#endif /* ALLOW_NONHYDROSTATIC */ + +c ELSE +C---- use constant viscosity (useVariableVisc=F): +c DO j=1-OLy,sNy+OLy +c DO i=1-OLx,sNx+OLx +c viscAh_D(i,j) = viscAhD +c viscAh_Z(i,j) = viscAhZ +c viscA4_D(i,j) = viscA4D +c viscA4_Z(i,j) = viscA4Z +c ENDDO +c ENDDO +C---- variable/constant viscosity : end if/else block +c ENDIF + +#ifdef ALLOW_DIAGNOSTICS + IF (useDiagnostics) THEN + CALL DIAGNOSTICS_FILL(viscAh_D,'VISCAHD ',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid) + + CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscA4_DMax,'VA4DMAX ',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscAh_ZMax,'VAHZMAX ',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscA4_ZMax,'VA4ZMAX ',k,1,2,bi,bj,myThid) + + CALL DIAGNOSTICS_FILL(viscAh_DMin,'VAHDMIN ',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscA4_DMin,'VA4DMIN ',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscAh_ZMin,'VAHZMIN ',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscA4_ZMin,'VA4ZMIN ',k,1,2,bi,bj,myThid) + + CALL DIAGNOSTICS_FILL(viscAh_DLth,'VAHDLTH ',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscA4_DLth,'VA4DLTH ',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscAh_ZLth,'VAHZLTH ',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscA4_ZLth,'VA4ZLTH ',k,1,2,bi,bj,myThid) + + CALL DIAGNOSTICS_FILL(viscAh_DLthD,'VAHDLTHD', + & k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscA4_DLthD,'VA4DLTHD', + & k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscAh_ZLthD,'VAHZLTHD', + & k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscA4_ZLthD,'VA4ZLTHD', + & k,1,2,bi,bj,myThid) +#ifdef ALLOW_LEITH_QG + CALL DIAGNOSTICS_FILL(viscAh_DLthQG,'VAHDLTHQ', + & k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscAh_ZLthQG,'VAHZLTHQ', + & k,1,2,bi,bj,myThid) +#endif + CALL DIAGNOSTICS_FILL(viscAh_DSmg,'VAHDSMAG',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscA4_DSmg,'VA4DSMAG',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscAh_ZSmg,'VAHZSMAG',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(viscA4_ZSmg,'VA4ZSMAG',k,1,2,bi,bj,myThid) + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ + + RETURN + END diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/packages.conf b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/packages.conf new file mode 100755 index 0000000..f3247d4 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/packages.conf @@ -0,0 +1,7 @@ +gfd +exch2 +diagnostics +ggl90 +gmredi +salt_plume +seaice diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_advdiff.F b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_advdiff.F new file mode 100755 index 0000000..de41a95 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_advdiff.F @@ -0,0 +1,771 @@ +#include "SEAICE_OPTIONS.h" +#ifdef ALLOW_AUTODIFF +# include "AUTODIFF_OPTIONS.h" +#endif + +CBOP +C !ROUTINE: SEAICE_ADVDIFF + +C !INTERFACE: ========================================================== + SUBROUTINE SEAICE_ADVDIFF( + U uc, vc, + I myTime, myIter, myThid ) + +C !DESCRIPTION: \bv +C *===========================================================* +C | SUBROUTINE SEAICE_ADVDIFF +C | o driver for different advection routines +C | calls an adaption of gad_advection to call different +C | advection routines of pkg/generic_advdiff +C *===========================================================* +C \ev + +C !USES: =============================================================== + IMPLICIT NONE + +C === Global variables === +#include "SIZE.h" +#include "EEPARAMS.h" +#include "PARAMS.h" +#include "GRID.h" +#include "SEAICE_SIZE.h" +#include "SEAICE_PARAMS.h" +#include "SEAICE.h" +#include "SEAICE_TRACER.h" +#ifdef HACK_FOR_GMAO_CPL +# include "SEAICE_LAYERS.h" +#endif +#ifdef ALLOW_AUTODIFF_TAMC +# include "tamc.h" +#endif + +C !INPUT/OUTPUT PARAMETERS: =================================================== +C === Routine arguments === +C uc/vc :: current ice velocity on C-grid; +C :: C-Grid : Input only ; B-grid : Output only +C myTime :: current time in simulation +C myIter :: current iteration number in simulation +C myThid :: my Thread Id number + _RL uc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL vc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL myTime + INTEGER myIter + INTEGER myThid + +C !LOCAL VARIABLES: ==================================================== +C === Local variables === +C i,j,bi,bj :: Loop counters +C it :: Loop counter for ice thickness categories +C uTrans :: volume transport, x direction +C vTrans :: volume transport, y direction +C afx :: horizontal advective flux, x direction +C afy :: horizontal advective flux, y direction +C gFld :: tendency of seaice field +C xA,yA :: "areas" of X and Y face of tracer cells + INTEGER i, j, bi, bj +#ifdef SEAICE_ITD + INTEGER it +#endif /* SEAICE_ITD */ +#ifdef HACK_FOR_GMAO_CPL + INTEGER l, n +#endif /* HACK_FOR_GMAO_CPL */ +#ifdef ALLOW_AUTODIFF_TAMC + INTEGER itmpkey +#endif /* ALLOW_AUTODIFF_TAMC */ +#ifdef ALLOW_SITRACER + _RL hEffNm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL areaNm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + INTEGER iTr, SEAICEadvSchSItr + _RL SEAICEdiffKhSItr + _RL SItrExt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL tmpscal1, tmpscal2 +# ifdef ALLOW_SITRACER_ADVCAP + _RL SItrPrev (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) +# endif +# ifdef ALLOW_SITRACER_DEBUG_DIAG + _RL DIAGarray (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) +# endif +#endif /* ALLOW_SITRACER */ + _RL fldNm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL afx (1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL afy (1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL gFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL recip_heff(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + CHARACTER*(MAX_LEN_MBUF) msgBuf +CEOP +#ifdef HACK_FOR_GMAO_CPL + INTEGER SEAICEadvSchQice, SEAICEadvSchQsnow + INTEGER SEAICEadvSchMltPd + SEAICEadvSchQice = 0 + SEAICEadvSchQsnow = 0 + SEAICEadvSchMltPd = 0 + SEAICEadvSchQice = SEAICEadvSchHeff + SEAICEadvSchQsnow = SEAICEadvSchSnow +c SEAICEadvSchMltPd = SEAICEadvSchSnow +#endif /* HACK_FOR_GMAO_CPL */ + +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| + +C-- make a local copy of the velocities for compatibility with B-grid +C-- alternatively interpolate to C-points if necessary + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) +#ifndef SEAICE_CGRID /* not SEAICE_CGRID = BGRID */ +C-- hack to ensure backward compatibility: +C average B-grid seaice velocity to C-grid + DO j=1-OLy,sNy+OLy-1 + DO i=1-OLx,sNx+OLx-1 + uc(i,j,bi,bj)=.5 _d 0*(UICE(i,j,bi,bj)+UICE(i,j+1,bi,bj)) + vc(i,j,bi,bj)=.5 _d 0*(VICE(i,j,bi,bj)+VICE(i+1,j,bi,bj)) + ENDDO + ENDDO +#endif /* SEAICE_CGRID */ +C- compute cell areas used by all tracers + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + xA(i,j,bi,bj) = _dyG(i,j,bi,bj)*SIMaskU(i,j,bi,bj) + yA(i,j,bi,bj) = _dxG(i,j,bi,bj)*SIMaskV(i,j,bi,bj) + ENDDO + ENDDO + ENDDO + ENDDO + +#ifndef SEAICE_CGRID +C Do we need this? I am afraid so. + CALL EXCH_UV_XY_RL(uc,vc,.TRUE.,myThid) +#endif /* not SEAICE_CGRID */ + +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE uc = comlev1, key = ikey_dynamics, kind=isbyte +CADJ STORE vc = comlev1, key = ikey_dynamics, kind=isbyte +#endif /* ALLOW_AUTODIFF_TAMC */ + + IF ( SEAICEmultiDimAdvection ) THEN +#ifdef ALLOW_GENERIC_ADVDIFF + + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) +C--- loops on tile indices bi,bj + +#ifdef ALLOW_AUTODIFF_TAMC +C Initialise for TAF + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + gFld(i,j) = 0. _d 0 + ENDDO + ENDDO +C + act1 = bi - myBxLo(myThid) + max1 = myBxHi(myThid) - myBxLo(myThid) + 1 + act2 = bj - myByLo(myThid) + max2 = myByHi(myThid) - myByLo(myThid) + 1 + act3 = myThid - 1 + max3 = nTx*nTy + act4 = ikey_dynamics - 1 + itmpkey = (act1 + 1) + act2*max1 + & + act3*max1*max2 + & + act4*max1*max2*max3 +CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key=itmpkey, kind=isbyte +CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key=itmpkey, kind=isbyte +CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key=itmpkey, kind=isbyte +# ifdef SEAICE_VARIABLE_SALINITY +CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj, key=itmpkey, kind=isbyte +# endif +#endif /* ALLOW_AUTODIFF_TAMC */ + + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx +#ifdef ALLOW_SITRACER + hEffNm1(i,j,bi,bj) = HEFF(i,j,bi,bj) + areaNm1(i,j,bi,bj) = AREA(i,j,bi,bj) +#endif + recip_heff(i,j) = 1. _d 0 + ENDDO + ENDDO + +C- Calculate "volume transports" through tracer cell faces. + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + uTrans(i,j) = uc(i,j,bi,bj)*xA(i,j,bi,bj) + vTrans(i,j) = vc(i,j,bi,bj)*yA(i,j,bi,bj) + ENDDO + ENDDO + +#ifdef SEAICE_ITD +C-- Effective Thickness (Volume) + IF ( SEAICEadvHeff ) THEN + DO it=1,SEAICE_multDim + CALL SEAICE_ADVECTION( + I GAD_HEFF, SEAICEadvSchHeff, + I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj), + I uTrans, vTrans, HEFFITD(1-OLx,1-OLy,it,bi,bj), + I recip_heff, + O gFld, afx, afy, + I bi, bj, myTime, myIter, myThid ) +C- Add tendency due to diffusion + IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) + & CALL SEAICE_DIFFUSION( + I GAD_HEFF, SEAICEdiffKhHeff, ONE, + I HEFFITD(1-OLx,1-OLy,it,bi,bj), HEFFM, + I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj), + U gFld, + I bi, bj, myTime, myIter, myThid ) +C now do the "explicit" time step + DO j=1,sNy + DO i=1,sNx + HEFFITD(i,j,it,bi,bj) = HEFFM(i,j,bi,bj) * ( + & HEFFITD(i,j,it,bi,bj) + SEAICE_deltaTtherm * gFld(i,j) + & ) + ENDDO + ENDDO + ENDDO + ENDIF + +C-- Fractional area + IF ( SEAICEadvArea ) THEN + DO it=1,SEAICE_multDim + CALL SEAICE_ADVECTION( + I GAD_AREA, SEAICEadvSchArea, + I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj), + I uTrans, vTrans, AREAITD(1-OLx,1-OLy,it,bi,bj), + I recip_heff, + O gFld, afx, afy, + I bi, bj, myTime, myIter, myThid ) +C- Add tendency due to diffusion + IF ( SEAICEdiffKhArea .GT. 0. _d 0 ) + & CALL SEAICE_DIFFUSION( + I GAD_AREA, SEAICEdiffKhArea, ONE, + I AREAITD(1-OLx,1-OLy,it,bi,bj), HEFFM, + I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj), + U gFld, + I bi, bj, myTime, myIter, myThid ) +C now do the "explicit" time step + DO j=1,sNy + DO i=1,sNx + AREAITD(i,j,it,bi,bj) = HEFFM(i,j,bi,bj) * ( + & AREAITD(i,j,it,bi,bj) + SEAICE_deltaTtherm * gFld(i,j) + & ) + ENDDO + ENDDO + ENDDO +C open water fraction needs to be advected for the ridging scheme + CALL SEAICE_ADVECTION( + I GAD_AREA, SEAICEadvSchArea, + I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj), + I uTrans, vTrans, opnWtrFrac(1-OLx,1-OLy,bi,bj), recip_heff, + O gFld, afx, afy, + I bi, bj, myTime, myIter, myThid ) +C-- Add tendency due to diffusion + IF ( SEAICEdiffKhArea .GT. 0. _d 0 ) + & CALL SEAICE_DIFFUSION( + I GAD_AREA, SEAICEdiffKhArea, ONE, + I opnWtrFrac(1-OLx,1-OLy,bi,bj), HEFFM, + I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj), + U gFld, + I bi, bj, myTime, myIter, myThid ) +C now do the "explicit" time step + DO j=1,sNy + DO i=1,sNx + opnWtrFrac(i,j,bi,bj) = HEFFM(i,j,bi,bj) * ( + & opnWtrFrac(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j) + & ) + ENDDO + ENDDO + ENDIF + +C-- Effective Snow Thickness (Volume) + IF ( SEAICEadvSnow ) THEN + DO it=1,SEAICE_multDim + CALL SEAICE_ADVECTION( + I GAD_SNOW, SEAICEadvSchSnow, + I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj), + I uTrans, vTrans, HSNOWITD(1-OLx,1-OLy,it,bi,bj), + I recip_heff, + O gFld, afx, afy, + I bi, bj, myTime, myIter, myThid ) +C-- Add tendency due to diffusion + IF ( SEAICEdiffKhSnow .GT. 0. _d 0 ) + & CALL SEAICE_DIFFUSION( + I GAD_SNOW, SEAICEdiffKhSnow, ONE, + I HSNOWITD(1-OLx,1-OLy,it,bi,bj), HEFFM, + I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj), + U gFld, + I bi, bj, myTime, myIter, myThid ) +C now do the "explicit" time step + DO j=1,sNy + DO i=1,sNx + HSNOWITD(i,j,it,bi,bj) = HEFFM(i,j,bi,bj) * ( + & HSNOWITD(i,j,it,bi,bj) + SEAICE_deltaTtherm*gFld(i,j) + & ) + ENDDO + ENDDO + ENDDO + ENDIF + +#ifdef HACK_FOR_GMAO_CPL + IF ( SEAICEadvSchQice.NE.0 ) THEN + DO n=1,SEAICE_multDim + DO l=1,nIceLayers + CALL SEAICE_ADVECTION( + I GAD_QICE1, SEAICEadvSchQice, + I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj), uTrans, + I vTrans, SIqIce(1-OLx,1-OLy,l,n,bi,bj), recip_heff, + O gFld, afx, afy, + I bi, bj, myTime, myIter, myThid ) + IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) THEN +C- Add tendency due to diffusion + CALL SEAICE_DIFFUSION( + I GAD_QICE1, SEAICEdiffKhHeff, oneRL, + I SIqIce(1-OLx,1-OLy,l,n,bi,bj), HEFFM, + I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj), + U gFld, + I bi, bj, myTime, myIter, myThid ) + ENDIF +C now do the "explicit" time step + DO j=1,sNy + DO i=1,sNx + SIqIce(i,j,l,n,bi,bj) = ( SIqIce(i,j,l,n,bi,bj) + & + SEAICE_deltaTtherm * gFld(i,j) + & ) * HEFFM(i,j,bi,bj) + ENDDO + ENDDO + ENDDO + ENDDO + ENDIF +C-- + IF ( SEAICEadvSchQsnow.NE.0 ) THEN + DO n=1,SEAICE_multDim + DO l=1,nSnowLayers + CALL SEAICE_ADVECTION( + I GAD_QICE2, SEAICEadvSchQsnow, + I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj), uTrans, + I vTrans, SIqSnow(1-OLx,1-OLy,l,n,bi,bj), recip_heff, + O gFld, afx, afy, + I bi, bj, myTime, myIter, myThid ) + IF ( SEAICEdiffKhSnow .GT. 0. _d 0 ) THEN +C- Add tendency due to diffusion + CALL SEAICE_DIFFUSION( + I GAD_QICE2, SEAICEdiffKhSnow, oneRL, + I SIqSnow(1-OLx,1-OLy,l,n,bi,bj), HEFFM, + I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj), + U gFld, + I bi, bj, myTime, myIter, myThid ) + ENDIF +C now do the "explicit" time step + DO j=1,sNy + DO i=1,sNx + SIqSnow(i,j,l,n,bi,bj) = ( SIqSnow(i,j,l,n,bi,bj) + & + SEAICE_deltaTtherm * gFld(i,j) + & ) * HEFFM(i,j,bi,bj) + ENDDO + ENDDO + ENDDO + ENDDO + ENDIF +C-- + IF ( SEAICEadvSchMltPd.NE.0 ) THEN + DO n=1,SEAICE_multDim + CALL SEAICE_ADVECTION( + I GAD_QICE2, SEAICEadvSchMltPd, + I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj), uTrans, + I vTrans, SImeltPd(1-OLx,1-OLy,n,bi,bj), recip_heff, + O gFld, afx, afy, + I bi, bj, myTime, myIter, myThid ) + IF ( SEAICEdiffKhSnow .GT. 0. _d 0 ) THEN +C- Add tendency due to diffusion + CALL SEAICE_DIFFUSION( + I GAD_QICE2, SEAICEdiffKhSnow, oneRL, + I SImeltPd(1-OLx,1-OLy,n,bi,bj), HEFFM, + I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj), + U gFld, + I bi, bj, myTime, myIter, myThid ) + ENDIF +C now do the "explicit" time step + DO j=1,sNy + DO i=1,sNx + SImeltPd(i,j,n,bi,bj) = ( SImeltPd(i,j,n,bi,bj) + & + SEAICE_deltaTtherm * gFld(i,j) + & ) * HEFFM(i,j,bi,bj) + ENDDO + ENDDO + ENDDO + ENDIF +#endif /* HACK_FOR_GMAO_CPL */ + +C update mean ice thickness HEFF and total ice concentration AREA +C to match single category values +C (necessary here because updated HEFF is used below for SItracer) + CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid) + +#else /* not SEAICE_ITD */ +C-- Effective Thickness (Volume) + IF ( SEAICEadvHeff ) THEN + CALL SEAICE_ADVECTION( + I GAD_HEFF, SEAICEadvSchHeff, + I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj), + I uTrans, vTrans, HEFF(1-OLx,1-OLy,bi,bj), recip_heff, + O gFld, afx, afy, + I bi, bj, myTime, myIter, myThid ) + IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) THEN +C- Add tendency due to diffusion + CALL SEAICE_DIFFUSION( + I GAD_HEFF, SEAICEdiffKhHeff, ONE, + I HEFF(1-OLx,1-OLy,bi,bj), HEFFM, + I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj), + U gFld, + I bi, bj, myTime, myIter, myThid ) + ENDIF +C now do the "explicit" time step + DO j=1,sNy + DO i=1,sNx + HEFF(i,j,bi,bj) = HEFFM(i,j,bi,bj) * ( + & HEFF(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j) + & ) + ENDDO + ENDDO + ENDIF + +C-- Fractional area + IF ( SEAICEadvArea ) THEN + CALL SEAICE_ADVECTION( + I GAD_AREA, SEAICEadvSchArea, + I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj), + I uTrans, vTrans, AREA(1-OLx,1-OLy,bi,bj), recip_heff, + O gFld, afx, afy, + I bi, bj, myTime, myIter, myThid ) + IF ( SEAICEdiffKhArea .GT. 0. _d 0 ) THEN +C- Add tendency due to diffusion + CALL SEAICE_DIFFUSION( + I GAD_AREA, SEAICEdiffKhArea, ONE, + I AREA(1-OLx,1-OLy,bi,bj), HEFFM, + I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj), + U gFld, + I bi, bj, myTime, myIter, myThid ) + ENDIF +C now do the "explicit" time step + DO j=1,sNy + DO i=1,sNx + AREA(i,j,bi,bj) = HEFFM(i,j,bi,bj) * ( + & AREA(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j) + & ) + ENDDO + ENDDO + ENDIF + +C-- Effective Snow Thickness (Volume) + IF ( SEAICEadvSnow ) THEN + CALL SEAICE_ADVECTION( + I GAD_SNOW, SEAICEadvSchSnow, + I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj), + I uTrans, vTrans, HSNOW(1-OLx,1-OLy,bi,bj), recip_heff, + O gFld, afx, afy, + I bi, bj, myTime, myIter, myThid ) + IF ( SEAICEdiffKhSnow .GT. 0. _d 0 ) THEN +C-- Add tendency due to diffusion + CALL SEAICE_DIFFUSION( + I GAD_SNOW, SEAICEdiffKhSnow, ONE, + I HSNOW(1-OLx,1-OLy,bi,bj), HEFFM, + I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj), + U gFld, + I bi, bj, myTime, myIter, myThid ) + ENDIF +C now do the "explicit" time step + DO j=1,sNy + DO i=1,sNx + HSNOW(i,j,bi,bj) = HEFFM(i,j,bi,bj) * ( + & HSNOW(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j) + & ) + ENDDO + ENDDO + ENDIF +#endif /* SEAICE_ITD */ + +#ifdef SEAICE_VARIABLE_SALINITY +C-- Effective Sea Ice Salinity (Mass of salt) + IF ( SEAICEadvSalt ) THEN + CALL SEAICE_ADVECTION( + I GAD_SALT, SEAICEadvSchSalt, + I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj), + I uTrans, vTrans, HSALT(1-OLx,1-OLy,bi,bj), recip_heff, + O gFld, afx, afy, + I bi, bj, myTime, myIter, myThid ) + IF ( SEAICEdiffKhSalt .GT. 0. _d 0 ) THEN +C-- Add tendency due to diffusion + CALL SEAICE_DIFFUSION( + I GAD_SALT, SEAICEdiffKhSalt, ONE, + I HSALT(1-OLx,1-OLy,bi,bj), HEFFM, + I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj), + U gFld, + I bi, bj, myTime, myIter, myThid ) + ENDIF +C now do the "explicit" time step + DO j=1,sNy + DO i=1,sNx + HSALT(i,j,bi,bj) = HEFFM(i,j,bi,bj) * ( + & HSALT(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j) + & ) + ENDDO + ENDDO + ENDIF +#endif /* SEAICE_VARIABLE_SALINITY */ + +#ifdef ALLOW_SITRACER +C-- Sea Ice Tracers + DO iTr = 1, SItrNumInUse + IF ( (SEAICEadvHEFF.AND.(SItrMate(iTr).EQ.'HEFF')).OR. + & (SEAICEadvAREA.AND.(SItrMate(iTr).EQ.'AREA')) ) THEN +C-- scale to effective value + IF (SItrMate(iTr).EQ.'HEFF') THEN + SEAICEadvSchSItr=SEAICEadvSchHEFF + SEAICEdiffKhSItr=SEAICEdiffKhHEFF + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + SItrExt(i,j,bi,bj) = HEFFM(i,j,bi,bj) * + & SItracer(i,j,bi,bj,iTr) * hEffNm1(i,j,bi,bj) + ENDDO + ENDDO +c TAF? ELSEIF (SItrMate(iTr).EQ.'AREA') THEN + ELSE + SEAICEadvSchSItr=SEAICEadvSchAREA + SEAICEdiffKhSItr=SEAICEdiffKhAREA + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + SItrExt(i,j,bi,bj) = HEFFM(i,j,bi,bj) * + & SItracer(i,j,bi,bj,iTr) * areaNm1(i,j,bi,bj) + ENDDO + ENDDO + ENDIF +C-- store a couple things + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx +#ifdef ALLOW_SITRACER_ADVCAP +C-- store previous value for spurious maxima treament + SItrPrev(i,j,bi,bj)=SItracer(i,j,bi,bj,iTr) +#endif +#ifdef ALLOW_SITRACER_DEBUG_DIAG + diagArray(I,J,2+(iTr-1)*5) = SItrExt(i,j,bi,bj) +#endif + ENDDO + ENDDO +C-- compute advective tendency + CALL SEAICE_ADVECTION( + I GAD_SITR+iTr-1, SEAICEadvSchSItr, + I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj), + I uTrans, vTrans, SItrExt(1-OLx,1-OLy,bi,bj), + I recip_heff, + O gFld, afx, afy, + I bi, bj, myTime, myIter, myThid ) + IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) THEN +C-- add diffusive tendency + CALL SEAICE_DIFFUSION( + I GAD_SITR+iTr-1, SEAICEdiffKhSItr, ONE, + I SItrExt(1-OLx,1-OLy,bi,bj), HEFFM, + I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj), + U gFld, + I bi, bj, myTime, myIter, myThid ) + ENDIF +C-- apply tendency + DO j=1,sNy + DO i=1,sNx + SItrExt(i,j,bi,bj) = HEFFM(i,j,bi,bj) * ( + & SItrExt(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j) ) + ENDDO + ENDDO +C-- scale back to actual value, or move effective value to ocean bucket + IF (SItrMate(iTr).EQ.'HEFF') THEN + DO j=1,sNy + DO i=1,sNx + if (HEFF(I,J,bi,bj).GE.siEps) then + SItracer(i,j,bi,bj,iTr)=SItrExt(i,j,bi,bj)/HEFF(I,J,bi,bj) + SItrBucket(i,j,bi,bj,iTr)=0. _d 0 + else + SItracer(i,j,bi,bj,iTr)=0. _d 0 + SItrBucket(i,j,bi,bj,iTr)=SItrExt(i,j,bi,bj) + endif +#ifdef ALLOW_SITRACER_ADVCAP +C hack to try avoid 'spontaneous generation' of maxima, which supposedly would +C occur less frequently if we advected SItr with uXheff instead SItrXheff with u + tmpscal1=max(SItrPrev(i,j,bi,bj), + & SItrPrev(i+1,j,bi,bj),SItrPrev(i-1,j,bi,bj), + & SItrPrev(i,j+1,bi,bj),SItrPrev(i,j-1,bi,bj)) + tmpscal2=MAX(ZERO,SItracer(i,j,bi,bj,iTr)-tmpscal1) + SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal2 + SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr) + & +tmpscal2*HEFF(I,J,bi,bj) +#endif +C treat case of potential negative value + if (HEFF(I,J,bi,bj).GE.siEps) then + tmpscal1=MIN(0. _d 0,SItracer(i,j,bi,bj,iTr)) + SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal1 + SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr) + & +HEFF(I,J,bi,bj)*tmpscal1 + endif +#ifdef ALLOW_SITRACER_DEBUG_DIAG + diagArray(I,J,1+(iTr-1)*5)= - SItrBucket(i,j,bi,bj,iTr) + & *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce + tmpscal1= ( HEFF(I,J,bi,bj)*SItracer(i,j,bi,bj,iTr) + & + SItrBucket(i,j,bi,bj,iTr) )*HEFFM(I,J,bi,bj) + diagArray(I,J,2+(iTr-1)*5)= tmpscal1-diagArray(I,J,2+(iTr-1)*5) + diagArray(I,J,3+(iTr-1)*5)=HEFFM(i,j,bi,bj) * + & SEAICE_deltaTtherm * gFld(i,j) +#endif + ENDDO + ENDDO +c TAF? ELSEIF (SItrMate(iTr).EQ.'AREA') THEN + ELSE + DO j=1,sNy + DO i=1,sNx + if (AREA(I,J,bi,bj).GE.SEAICE_area_floor) then + SItracer(i,j,bi,bj,iTr)=SItrExt(i,j,bi,bj)/AREA(I,J,bi,bj) + else + SItracer(i,j,bi,bj,iTr)=0. _d 0 + endif + SItrBucket(i,j,bi,bj,iTr)=0. _d 0 +#ifdef ALLOW_SITRACER_ADVCAP + tmpscal1=max(SItrPrev(i,j,bi,bj), + & SItrPrev(i+1,j,bi,bj),SItrPrev(i-1,j,bi,bj), + & SItrPrev(i,j+1,bi,bj),SItrPrev(i,j-1,bi,bj)) + tmpscal2=MAX(ZERO,SItracer(i,j,bi,bj,iTr)-tmpscal1) + SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal2 +#endif +C treat case of potential negative value + if (AREA(I,J,bi,bj).GE.SEAICE_area_floor) then + tmpscal1=MIN(0. _d 0,SItracer(i,j,bi,bj,iTr)) + SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal1 + endif +#ifdef ALLOW_SITRACER_DEBUG_DIAG + diagArray(I,J,1+(iTr-1)*5)= 0. _d 0 + diagArray(I,J,2+(iTr-1)*5)= - diagArray(I,J,2+(iTr-1)*5) + & + AREA(I,J,bi,bj)*SItracer(i,j,bi,bj,iTr)*HEFFM(I,J,bi,bj) + diagArray(I,J,3+(iTr-1)*5)=HEFFM(i,j,bi,bj) * + & SEAICE_deltaTtherm * gFld(i,j) +#endif + ENDDO + ENDDO + ENDIF +C-- + ENDIF + ENDDO +#ifdef ALLOW_SITRACER_DEBUG_DIAG +c CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG2 ',0,Nr,2,bi,bj,myThid) +#endif +#endif /* ALLOW_SITRACER */ + +C--- end bi,bj loops + ENDDO + ENDDO + +#else /* not ALLOW_GENERIC_ADVDIFF */ + WRITE(msgBuf,'(2A)') + & 'SEAICE_ADVDIFF: cannot use SEAICEmultiDimAdvection', + & ' without pkg/generic_advdiff' + CALL PRINT_ERROR( msgBuf , myThid ) + WRITE(msgBuf,'(2A)') 'SEAICE_ADVDIFF: ', + & 'Re-compile with pkg "generic_advdiff" in packages.conf' + CALL PRINT_ERROR( msgBuf , myThid ) + CALL ALL_PROC_DIE( myThid ) + STOP 'ABNORMAL END: S/R SEAICE_ADVDIFF' +#endif /* ALLOW_GENERIC_ADVDIFF */ + ELSE +C-- if not multiDimAdvection +#ifdef SEAICE_ITD +C just for safety + WRITE(msgBuf,'(2A)') 'SEAICE_ADVDIFF: ', + & 'ITD with SEAICEmultiDimAdvection=.False. is not allowed,' + CALL PRINT_ERROR( msgBuf , myThid ) + WRITE(msgBuf,'(2A)') 'SEAICE_ADVDIFF: ', + & 'use a multidimensional advection scheme (in data.seaice)' + CALL PRINT_ERROR( msgBuf , myThid ) + CALL ALL_PROC_DIE( myThid ) + STOP 'ABNORMAL END: S/R SEAICE_ADVDIFF' +#endif /* SEAICE_ITD */ + + IF ( SEAICEadvHEff ) THEN +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE heff = comlev1, key = ikey_dynamics, kind=isbyte +#endif + CALL ADVECT( uc, vc, hEff, fldNm1, HEFFM, myThid ) + IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) THEN +C- Add tendency due to diffusion + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + CALL SEAICE_DIFFUSION( + I GAD_HEFF, SEAICEdiffKhHeff, SEAICE_deltaTtherm, + I fldNm1(1-OLx,1-OLy,bi,bj), HEFFM, + I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj), + U HEFF(1-OLx,1-OLy,bi,bj), + I bi, bj, myTime, myIter, myThid ) + ENDDO + ENDDO + ENDIF + ENDIF + IF ( SEAICEadvArea ) THEN +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE area = comlev1, key = ikey_dynamics, kind=isbyte +#endif + CALL ADVECT( uc, vc, area, fldNm1, HEFFM, myThid ) + IF ( SEAICEdiffKhArea .GT. 0. _d 0 ) THEN +C- Add tendency due to diffusion + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + CALL SEAICE_DIFFUSION( + I GAD_AREA, SEAICEdiffKhArea, SEAICE_deltaTtherm, + I fldNm1(1-OLx,1-OLy,bi,bj), HEFFM, + I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj), + U Area(1-OLx,1-OLy,bi,bj), + I bi, bj, myTime, myIter, myThid ) + ENDDO + ENDDO + ENDIF + ENDIF + IF ( SEAICEadvSnow ) THEN +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE hsnow = comlev1, key = ikey_dynamics, kind=isbyte +#endif + CALL ADVECT( uc, vc, HSNOW, fldNm1, HEFFM, myThid ) + IF ( SEAICEdiffKhSnow .GT. 0. _d 0 ) THEN +C- Add tendency due to diffusion + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + CALL SEAICE_DIFFUSION( + I GAD_SNOW, SEAICEdiffKhSnow, SEAICE_deltaTtherm, + I fldNm1(1-OLx,1-OLy,bi,bj), HEFFM, + I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj), + U HSNOW(1-OLx,1-OLy,bi,bj), + I bi, bj, myTime, myIter, myThid ) + ENDDO + ENDDO + ENDIF + ENDIF + +#ifdef SEAICE_VARIABLE_SALINITY + IF ( SEAICEadvSalt ) THEN +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE hsalt = comlev1, key = ikey_dynamics, kind=isbyte +#endif + CALL ADVECT( uc, vc, HSALT, fldNm1, HEFFM, myThid ) + IF ( SEAICEdiffKhSalt .GT. 0. _d 0 ) THEN +C- Add tendency due to diffusion + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + CALL SEAICE_DIFFUSION( + I GAD_SALT, SEAICEdiffKhSalt, SEAICE_deltaTtherm, + I fldNm1(1-OLx,1-OLy,bi,bj), HEFFM, + I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj), + U HSALT(1-OLx,1-OLy,bi,bj), + I bi, bj, myTime, myIter, myThid ) + ENDDO + ENDDO + ENDIF + ENDIF +#endif /* SEAICE_VARIABLE_SALINITY */ + +C-- end if multiDimAdvection + ENDIF + + RETURN + END diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_diag_init_add.h b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_diag_init_add.h new file mode 100755 index 0000000..453b74d --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_diag_init_add.h @@ -0,0 +1,120 @@ + + diagName = 'CPLoWGHT' + diagTitle = 'grid-cell Ocean fraction from GEOS' + diagUnits = '1 ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + +#ifdef SEAICE_ITD + + diagName = 'SItIces ' + diagTitle = 'Surface Temperature over Seaice for each category' + diagUnits = 'K ' + diagCode = 'SM MX ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + CALL DIAGNOSTICS_SETKLEV( diagName, nITD, myThid ) + + diagName = 'SIqIce ' + diagTitle = 'SEAICE enthalpy for each layer and category' + diagUnits = 'J/m^2 ' + diagCode = 'SM MX ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + numArea = nIceLayers*nITD + CALL DIAGNOSTICS_SETKLEV( diagName, numArea, myThid ) + + diagName = 'SIqSnow ' + diagTitle = 'Snow enthalpy for each layer and category' + diagUnits = 'J/m^2 ' + diagCode = 'SM MX ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + numArea = nSnowLayers*nITD + CALL DIAGNOSTICS_SETKLEV( diagName, numArea, myThid ) + + diagName = 'SImeltPd' + diagTitle = 'Melt Pond volume for each category' + diagUnits = 'm ' + diagCode = 'SM MX ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + CALL DIAGNOSTICS_SETKLEV( diagName, nITD, myThid ) + + diagName = 'SIiceAge' + diagTitle = 'Seaice Age for each category' + diagUnits = 's ' + diagCode = 'SM MX ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + CALL DIAGNOSTICS_SETKLEV( diagName, nITD, myThid ) + +C- Advection Increment: + diagName = 'SI_dArea' + diagTitle = 'Seaice fraction Advection Increment per cat.' + diagUnits = '1 ' + diagCode = 'SM MX ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + CALL DIAGNOSTICS_SETKLEV( diagName, nITD, myThid ) + + diagName = 'SI_dHeff' + diagTitle = 'Seaice thickness Advection Increment per cat.' + diagUnits = 'm ' + diagCode = 'SM MX ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + CALL DIAGNOSTICS_SETKLEV( diagName, nITD, myThid ) + + diagName = 'SI_dHsnw' + diagTitle = 'Snow thickness Advection Increment per cat.' + diagUnits = 'm ' + diagCode = 'SM MX ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + CALL DIAGNOSTICS_SETKLEV( diagName, nITD, myThid ) + + diagName = 'SI_dTIce' + diagTitle = 'Seaice surf. temp. Advection Increment per cat.' + diagUnits = 'K ' + diagCode = 'SM MX ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + CALL DIAGNOSTICS_SETKLEV( diagName, nITD, myThid ) + + diagName = 'SI_dQIce' + diagTitle = 'Seaice enthalpy Advect. Increment per layer and cat.' + diagUnits = 'J/m^2 ' + diagCode = 'SM MX ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + numArea = nIceLayers*nITD + CALL DIAGNOSTICS_SETKLEV( diagName, numArea, myThid ) + + diagName = 'SI_dQSnw' + diagTitle = 'Snow enthalpy Advect. Increment per layer and cat.' + diagUnits = 'J/m^2 ' + diagCode = 'SM MX ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + numArea = nSnowLayers*nITD + CALL DIAGNOSTICS_SETKLEV( diagName, numArea, myThid ) + + diagName = 'SI_dMPnd' + diagTitle = 'Melt Pond volume Advection Increment per cat.' + diagUnits = 'm ' + diagCode = 'SM MX ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + CALL DIAGNOSTICS_SETKLEV( diagName, nITD, myThid ) + + diagName = 'SI_dIcAg' + diagTitle = 'Seaice Age Advection Increment per cat.' + diagUnits = 's ' + diagCode = 'SM MX ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + CALL DIAGNOSTICS_SETKLEV( diagName, nITD, myThid ) + +#endif /* SEAICE_ITD */ diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_diagnostics_init.F b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_diagnostics_init.F new file mode 100755 index 0000000..d95482b --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_diagnostics_init.F @@ -0,0 +1,838 @@ +#include "SEAICE_OPTIONS.h" + +C-- File seaice_diagnostics_init.F: Routines initialize SEAICE diagnostics +C-- Contents +C-- o SEAICE_DIAGNOSTICS_INIT +C-- o SEAICE_DIAG_SUFX + +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| +CBOP +C !ROUTINE: SEAICE_DIAGNOSTICS_INIT +C !INTERFACE: + SUBROUTINE SEAICE_DIAGNOSTICS_INIT( myThid ) + +C !DESCRIPTION: \bv +C *==========================================================* +C | SUBROUTINE SEAICE_DIAGNOSTICS_INIT +C | o Routine to initialize list of all available diagnostics +C | for SEAICE package +C *==========================================================* +C \ev +C !USES: + IMPLICIT NONE + +C === Global variables === +#include "EEPARAMS.h" +#include "SIZE.h" +#include "SEAICE_SIZE.h" +#include "SEAICE_PARAMS.h" +#include "SEAICE_TRACER.h" + +C !INPUT/OUTPUT PARAMETERS: +C === Routine arguments === +C myThid :: my Thread Id number + INTEGER myThid +CEOP + +#ifdef ALLOW_DIAGNOSTICS +C !LOCAL VARIABLES: +C === Local variables === +C diagNum :: diagnostics number in the (long) list of available diag. +C diagMate :: diag. mate number in the (long) list of available diag. +C diagName :: local short name (8c) of a diagnostics +C diagCode :: local parser field with characteristics of the diagnostics +C see head of S/R DIAGNOSTICS_INIT_EARLY or DIAGNOSTICS_MAIN_INIT +C for a list of options +C diagUnits :: local string (16c): physical units of a diagnostic field +C diagTitle :: local string (80c): description of field in diagnostic + INTEGER diagNum + INTEGER diagMate + CHARACTER*8 diagName + CHARACTER*16 diagCode + CHARACTER*16 diagUnits + CHARACTER*(80) diagTitle + +#ifdef ALLOW_SITRACER + INTEGER iTr, ilnb, numMateTr, numMateTrPreTh + CHARACTER*8 locUnitTr + CHARACTER*30 locNameTr +#endif + INTEGER numArea,numAreaPreTh,numHeff,numHeffPreTh + CHARACTER*9 flxUnits + CHARACTER*15 locName + CHARACTER*4 SEAICE_DIAG_SUFX, diagSufx + EXTERNAL SEAICE_DIAG_SUFX +C Functions :: + INTEGER ILNBLNK + EXTERNAL ILNBLNK + +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| + +C=============== state variables ============ + + diagName = 'SIarea ' + diagTitle = 'SEAICE fractional ice-covered area [0 to 1]' + diagUnits = 'm^2/m^2 ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + numArea = diagNum + + diagName = 'SIareaPR' + diagTitle = 'SIarea preceeding ridging process' + diagUnits = 'm^2/m^2 ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIareaPT' + diagTitle = 'SIarea preceeding thermodynamic growth/melt' + diagUnits = 'm^2/m^2 ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + numAreaPreTh = diagNum + + diagName = 'SIheff ' + diagTitle = 'SEAICE effective ice thickness' + diagUnits = 'm ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + numHeff = diagNum + + diagName = 'SIheffPT' + diagTitle = 'SIheff preceeeding thermodynamic growth/melt' + diagUnits = 'm ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + numHeffPreTh = diagNum + + diagName = 'SIhsnow ' + diagTitle = 'SEAICE effective snow thickness' + diagUnits = 'm ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIhsnoPT' + diagTitle = 'SIhsnow preceeeding thermodynamic growth/melt' + diagUnits = 'm ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIhsalt ' + diagTitle = 'SEAICE effective salinity' + diagUnits = 'g/m^2 ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + +#ifdef ALLOW_SITRACER + DO iTr = 1, SItrNumInUse + +C-- Set default name & tracer Units: + WRITE(locNameTr,'(A,I4.4,A)') 'sea ice tracer no. ',iTr + if (SItrMate(iTr).EQ.'HEFF') then + locUnitTr = '(kg/kg) ' + numMateTr = numHeff + numMateTrPreTh = numHeffPreTh + else + locUnitTr = '(kg/m^2)' + numMateTr = numArea + numMateTrPreTh = numAreaPreTh + endif +C- use name & units from data.seaice : + ilnb = ILNBLNK(SItrUnit(iTr)) + IF ( ilnb.GE.1 ) THEN + ilnb = LEN(locUnitTr) + locUnitTr = SItrUnit(iTr)(1:ilnb) + ENDIF + ilnb = ILNBLNK(SItrNameLong(iTr)) + IF ( ilnb.GE.1 ) THEN + ilnb = MIN(LEN(locNameTr),ilnb) + WRITE(locNameTr,'(A)') SItrNameLong(iTr)(1:ilnb) + ELSE + ilnb = ILNBLNK(SItrName(iTr)) + IF ( ilnb.GE.1 ) THEN + ilnb = MIN(LEN(locNameTr),ilnb) + WRITE(locNameTr,'(2A)') SItrName(iTr)(1:ilnb),' tracer' + ENDIF + ENDIF + ilnb = MAX(ILNBLNK(locNameTr),1) + + WRITE(diagName,'(A4,I2.2,A2)') 'SItr',iTr,' ' + WRITE(diagTitle,'(4A)') locNameTr(1:ilnb), + & ' (associated with ',SItrMate(iTr),')' + diagUnits = locUnitTr//' ' + diagCode = 'SM C M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, diagName, + I diagCode, diagUnits, diagTitle, numMateTr, myThid ) + + WRITE(diagName,'(A4,I2.2,A2)') 'SItr',iTr,'PT' + WRITE(diagTitle,'(A4,I2.2,2A)') 'SItr',iTr, + & ' preceeeding thermodynamic growth/melt' + diagUnits = locUnitTr//' ' + diagCode = 'SM C M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, diagName, + I diagCode, diagUnits, diagTitle, numMateTrPreTh, myThid ) + + ENDDO +#endif + + diagName = 'SItices ' + diagTitle = 'Surface Temperature over Sea-Ice (area weighted)' + diagUnits = 'K ' + diagCode = 'SM C M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, numArea, myThid ) + + diagName = 'SIuice ' + diagTitle = 'SEAICE zonal ice velocity, >0 from West to East' + diagUnits = 'm/s ' +#ifdef SEAICE_CGRID + diagCode = 'UU M1 ' +#else + diagCode = 'UZ M1 ' +#endif + diagMate = diagNum + 2 + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + + diagName = 'SIvice ' + diagTitle = 'SEAICE merid. ice velocity, >0 from South to North' + diagUnits = 'm/s ' +#ifdef SEAICE_CGRID + diagCode = 'VV M1 ' +#else + diagCode = 'VZ M1 ' +#endif + diagMate = diagNum + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + +C=============== momentum, heat and fresh water forcing ============ + +C pkg/diagnostics SIfu and oceTAUX, dumpfreq FU, and tavefreq FUtave +C are identical but they differ from pkg/diagnostics EXFtaux, which +C is stress before impact of ice. Also when using exf bulk +C formulae, EXFtaux is defined on tracer rather than uvel points. +c diagName = 'SIfu ' +c diagTitle = 'SEAICE zonal surface wind stress, >0 increases uVel ' +c diagUnits = 'N/m^2 ' +c diagCode = 'UU U1 ' +c diagMate = diagNum + 2 +c CALL DIAGNOSTICS_ADDTOLIST( diagNum, +c I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + +C pkg/diagnostics SIfv and oceTAUY, dumpfreq FV, and tavefreq FVtave +C are identical but they differ from pkg/diagnostics EXFtauy, which +C is stress before impact of ice. Also when using exf bulk +C formulae, EXFtauy is defined on tracer rather than vvel points. +c diagName = 'SIfv ' +c diagTitle = 'SEAICE merid. surface wind stress, >0 increases vVel' +c diagUnits = 'N/m^2 ' +c diagCode = 'VV U1 ' +c diagMate = diagNum +c CALL DIAGNOSTICS_ADDTOLIST( diagNum, +c I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + + diagName = 'SItaux ' + diagTitle = 'SEAICE zonal surface wind stress, >0 increases uIce' + diagUnits = 'N/m^2 ' + diagCode = 'UU U1 ' + diagMate = diagNum + 2 + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + + diagName = 'SItauy ' + diagTitle = 'SEAICE merid surface wind stress, >0 increases vIce' + diagUnits = 'N/m^2 ' + diagCode = 'VV U1 ' + diagMate = diagNum + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + + diagName = 'SIatmTx ' + diagTitle = 'Zonal surface wind stress over Ocean+SeaIce' + diagUnits = 'N/m^2 ' + diagCode = 'UU U1 ' + diagMate = diagNum + 2 + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + + diagName = 'SIatmTy ' + diagTitle = 'Merid surface wind stress over Ocean+SeaIce' + diagUnits = 'N/m^2 ' + diagCode = 'VV U1 ' + diagMate = diagNum + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + +c diagName = 'SIuwind ' +c diagTitle = 'SEAICE zonal 10-m wind speed, >0 increases uVel' +c diagUnits = 'm/s ' +c diagCode = 'UM U1 ' +c diagMate = diagNum + 2 +c CALL DIAGNOSTICS_ADDTOLIST( diagNum, +c I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + +c diagName = 'SIvwind ' +c diagTitle = 'SEAICE meridional 10-m wind speed, >0 increases uVel' +c diagUnits = 'm/s ' +c diagCode = 'VM U1 ' +c diagMate = diagNum +c CALL DIAGNOSTICS_ADDTOLIST( diagNum, +c I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + +C SIqnet, Qnet, and QNETtave are identical. +C With #undef NONLIN_FRSURF SIqnet is identical to -(TFLUX-TRELAX). +C Except over land and under sea ice, SIqnet is also identical to +C EXFlwnet+EXFswnet-EXFhl-EXFhs. + diagName = 'SIqnet ' + diagTitle = 'Ocean surface heatflux, turb+rad, >0 decreases theta' + diagUnits = 'W/m^2 ' + diagCode = 'SM U1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + +C SIqsw, Qsw, and QSWtave are identical. +C Except under sea ice, SIqsw is also identical to EXFswnet. + diagName = 'SIqsw ' + diagTitle = 'Ocean surface shortwave radiat., >0 decreases theta' + diagUnits = 'W/m^2 ' + diagCode = 'SM U1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIatmQnt' + diagTitle = 'Net atmospheric heat flux, >0 decreases theta' + diagUnits = 'W/m^2 ' + diagCode = 'SM U1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SItflux ' + diagTitle = 'Same as TFLUX but incl seaice (>0 incr T decr H)' + diagUnits = 'W/m^2 ' + diagCode = 'SM U1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + +#ifndef SEAICE_DISABLE_HEATCONSFIX + diagName = 'SIaaflux' + diagTitle = 'conservative ocn<->seaice adv. heat flux adjust.' + diagUnits = 'W/m^2 ' + diagCode = 'SM U1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) +#endif + + diagName = 'SIhl ' + diagTitle = 'Latent heat flux into ocean, >0 increases theta' + diagUnits = 'W/m^2 ' + diagCode = 'SM U1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIqneto ' + diagTitle = 'Open Ocean Part of SIqnet, turb+rad, >0 decr theta' + diagUnits = 'W/m^2 ' + diagCode = 'SM U1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIqneti ' + diagTitle = 'Ice Covered Part of SIqnet, turb+rad, >0 decr theta' + diagUnits = 'W/m^2 ' + diagCode = 'SM U1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + +C pkg/diagnostics SIempmr, dumpfreq EmPmR, and tavefreq EmPmRtave +C are identical but they differ from pkg/diagnostics EXFempmr, which +C is EmPmR before impact of ice. + diagName = 'SIempmr ' + diagTitle = 'Ocean surface freshwater flux, > 0 increases salt' + diagUnits = 'kg/m^2/s ' + diagCode = 'SM U1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIatmFW ' + diagTitle = 'Net freshwater flux from atmosphere & land (+=down)' + diagUnits = 'kg/m^2/s ' + diagCode = 'SM U1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIsnPrcp' + diagTitle = 'Snow precip. (+=dw) over Sea-Ice (area weighted)' + diagUnits = 'kg/m^2/s ' + diagCode = 'SM U1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIfwSubl' + diagTitle ='Potential sublimation freshwater flux, >0 decr. ice' + diagUnits = 'kg/m^2/s ' + diagCode = 'SM U1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIacSubl' + diagTitle = 'Actual sublimation freshwater flux, >0 decr. ice' + diagUnits = 'kg/m^2/s ' + diagCode = 'SM U1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIrsSubl' + diagTitle = 'Residual subl. freshwater flux, >0 taken from ocn' + diagUnits = 'kg/m^2/s ' + diagCode = 'SM U1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIactLHF' + diagTitle = 'Actual latent heat flux over ice' + diagUnits = 'W/m^2 ' + diagCode = 'SM U1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SImaxLHF' + diagTitle = 'Maximum latent heat flux over ice' + diagUnits = 'W/m^2 ' + diagCode = 'SM U1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + +#ifdef ALLOW_SITRACER + DO iTr = 1, SItrNumInUse + IF (SItrMate(iTr).EQ.'HEFF') then +C-- Set default name & tracer Units: + WRITE(diagUnits,'(A)') 'kg/m^2/s' +C-- use units from data.seaice : + ilnb = ILNBLNK(SItrUnit(iTr)) + IF ( ilnb.GE.1 ) THEN + WRITE(diagUnits,'(2A)') SItrUnit(iTr)(1:ilnb),'.kg/m^2/s' + ENDIF +C-- + WRITE(diagName,'(A4,I2.2,A2)') 'SItr',iTr,'Fx' + WRITE(diagTitle,'(A4,I2.2,A)') 'SItr',iTr, + I ' flux out of ice pack (that may enter ocean)' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, diagName, + I diagCode, diagUnits, diagTitle, 0, myThid ) + + ENDIF + ENDDO +#endif + +C============== ice growth/melt ============== + + diagName = 'SIaQbOCN' + diagTitle = 'Potential HEFF rate of change by ocean ice flux' + diagUnits = 'm/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIaQbATC' + diagTitle = 'Potential HEFF rate of change by atm flux over ice' + diagUnits = 'm/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIaQbATO' + diagTitle = 'Potential HEFF rate of change by open ocn atm flux' + diagUnits = 'm/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIdHbOCN' + diagTitle = 'HEFF rate of change by ocean ice flux' + diagUnits = 'm/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIdSbATC' + diagTitle = 'HSNOW rate of change by atm flux over sea ice' + diagUnits = 'm/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIdSbOCN' + diagTitle = 'HSNOW rate of change by ocean ice flux' + diagUnits = 'm/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIdHbATC' + diagTitle = 'HEFF rate of change by atm flux over sea ice' + diagUnits = 'm/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIdHbATO' + diagTitle = 'HEFF rate of change by open ocn atm flux' + diagUnits = 'm/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIdHbFLO' + diagTitle = 'HEFF rate of change by flooding snow' + diagUnits = 'm/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) +#ifdef SEAICE_GREASE + + diagName = 'SIgrsLT ' + diagTitle = 'actual grease ice layer thickness' + diagUnits = 'm ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) +#endif + +C=============== expansion/contraction ============ + + diagName = 'SIdAbATO' + diagTitle = 'Potential AREA rate of change by open ocn atm flux' + diagUnits = 'm^2/m^2/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIdAbATC' + diagTitle = 'Potential AREA rate of change by atm flux over ice' + diagUnits = 'm^2/m^2/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIdAbOCN' + diagTitle = 'Potential AREA rate of change by ocean ice flux' + diagUnits = 'm^2/m^2/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIdA' + diagTitle = 'AREA rate of change (net)' + diagUnits = 'm^2/m^2/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + +C============== advection/diffusion ============ + +C effective thickness + flxUnits = '.m^2/s ' + locName = 'eff ice thickn ' + WRITE(diagUnits,'(2A)') 'm',flxUnits + diagSufx = SEAICE_DIAG_SUFX( GAD_HEFF, myThid ) + +C-- advective flux + diagName = 'ADVx'//diagSufx + diagTitle = 'Zonal Advective Flux of '//locName + diagCode = 'UU M1 ' + diagMate = diagNum + 2 + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + + diagName = 'ADVy'//diagSufx + diagTitle = 'Meridional Advective Flux of '//locName + diagCode = 'VV M1 ' + diagMate = diagNum + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + +C-- Diffusive flux: + diagName = 'DFxE'//diagSufx + diagTitle = 'Zonal Diffusive Flux of '//locName + diagCode = 'UU M1 ' + diagMate = diagNum + 2 + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + + diagName = 'DFyE'//diagSufx + diagTitle = 'Meridional Diffusive Flux of '//locName + diagCode = 'VV M1 ' + diagMate = diagNum + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + +C fractional ice covered area (ice concentration) + locName = 'fract area ' + WRITE(diagUnits,'(2A)') 'm^2/m^2',flxUnits + diagSufx = SEAICE_DIAG_SUFX( GAD_AREA, myThid ) + +C-- advective flux + diagName = 'ADVx'//diagSufx + diagTitle = 'Zonal Advective Flux of '//locName + diagCode = 'UU M1 ' + diagMate = diagNum + 2 + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + + diagName = 'ADVy'//diagSufx + diagTitle = 'Meridional Advective Flux of '//locName + diagCode = 'VV M1 ' + diagMate = diagNum + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + +C-- Diffusive flux: + diagName = 'DFxE'//diagSufx + diagTitle = 'Zonal Diffusive Flux of '//locName + diagCode = 'UU M1 ' + diagMate = diagNum + 2 + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + + diagName = 'DFyE'//diagSufx + diagTitle = 'Meridional Diffusive Flux of '//locName + diagCode = 'VV M1 ' + diagMate = diagNum + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + +C effective snow thickness + locName = 'eff snow thickn' + WRITE(diagUnits,'(2A)') 'm',flxUnits + diagSufx = SEAICE_DIAG_SUFX( GAD_SNOW, myThid ) + +C-- advective flux + diagName = 'ADVx'//diagSufx + diagTitle = 'Zonal Advective Flux of '//locName + diagCode = 'UU M1 ' + diagMate = diagNum + 2 + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + + diagName = 'ADVy'//diagSufx + diagTitle = 'Meridional Advective Flux of '//locName + diagCode = 'VV M1 ' + diagMate = diagNum + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + +C-- Diffusive flux: + diagName = 'DFxE'//diagSufx + diagTitle = 'Zonal Diffusive Flux of '//locName + diagCode = 'UU M1 ' + diagMate = diagNum + 2 + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + + diagName = 'DFyE'//diagSufx + diagTitle = 'Meridional Diffusive Flux of '//locName + diagCode = 'VV M1 ' + diagMate = diagNum + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + +C sea ice salinity + locName = 'seaice salinity' + WRITE(diagUnits,'(2A)') 'psu',flxUnits + diagSufx = SEAICE_DIAG_SUFX( GAD_SALT, myThid ) + +C-- advective flux + diagName = 'ADVx'//diagSufx + diagTitle = 'Zonal Advective Flux of '//locName + diagCode = 'UU M1 ' + diagMate = diagNum + 2 + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + + diagName = 'ADVy'//diagSufx + diagTitle = 'Meridional Advective Flux of '//locName + diagCode = 'VV M1 ' + diagMate = diagNum + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + +C-- Diffusive flux: + diagName = 'DFxE'//diagSufx + diagTitle = 'Zonal Diffusive Flux of '//locName + diagCode = 'UU M1 ' + diagMate = diagNum + 2 + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + + diagName = 'DFyE'//diagSufx + diagTitle = 'Meridional Diffusive Flux of '//locName + diagCode = 'VV M1 ' + diagMate = diagNum + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + +C-- effective thickness transport (centered in space, 1 time-step lag) + diagName = 'SIuheff ' + diagTitle = 'Zonal Transport of eff ice thickn (centered)' + diagUnits = 'm^2/s ' + diagCode = 'UU M1 ' + diagMate = diagNum + 2 + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + + diagName = 'SIvheff ' + diagTitle = 'Meridional Transport of eff ice thickn (centered)' + diagUnits = 'm^2/s ' + diagCode = 'VV M1 ' + diagMate = diagNum + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + +C=============== dynamics ============ + + diagName = 'SIpress ' + diagTitle = 'SEAICE strength (with upper and lower limit)' + diagUnits = 'N/m ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIzeta ' + diagTitle = 'SEAICE nonlinear bulk viscosity' + diagUnits = 'kg/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIeta ' + diagTitle = 'SEAICE nonlinear shear viscosity' + diagUnits = 'kg/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIsig1 ' + diagTitle = 'SEAICE normalized principle stress, component one' + diagUnits = 'no units ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIsig2 ' + diagTitle = 'SEAICE normalized principle stress, component two' + diagUnits = 'no units ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIshear ' + diagTitle = 'SEAICE shear deformation rate' + diagUnits = '1/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIdelta ' + diagTitle = 'SEAICE Delta deformation rate' + diagUnits = '1/s ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SItensil' + diagTitle = 'SEAICE maximal tensile strength' + diagUnits = 'N/m ' + diagCode = 'SM M1 ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + +#ifdef SEAICE_ITD +C=============== ice thickness categories ============ + + diagName = 'SIheffN ' + diagTitle = 'SEAICE effective ice thickness per category' + diagUnits = 'm ' + diagCode = 'SM MX ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + CALL DIAGNOSTICS_SETKLEV( diagName, nITD, myThid ) + + diagName = 'SIareaN ' + diagTitle = + I 'SEAICE fractional ice-covered area per category [0 to 1]' + diagUnits = 'm^2/m^2 ' + diagCode = 'SM MX ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + CALL DIAGNOSTICS_SETKLEV( diagName, nITD, myThid ) + + diagName = 'SIhsnowN' + diagTitle = 'SEAICE effective snow thickness per category' + diagUnits = 'm ' + diagCode = 'SM MX ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + CALL DIAGNOSTICS_SETKLEV( diagName, nITD, myThid ) +#endif + +#ifdef HACK_FOR_GMAO_CPL +# include "seaice_diag_init_add.h" +#endif + +#endif /* ALLOW_DIAGNOSTICS */ + + RETURN + END + +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| +CBOP 0 +C !ROUTINE: SEAICE_DIAG_SUFX + +C !INTERFACE: + CHARACTER*4 FUNCTION SEAICE_DIAG_SUFX( tracerId, myThid ) + +C !DESCRIPTION: +C *==========================================================* +C | FUNCTION SEAICE_DIAG_SUFX +C | o Return diagnostic suffix (4 character long) for the +C | "tracerId" tracer (used to build diagnostic names). +C *==========================================================* + +C !USES: + IMPLICIT NONE +#include "EEPARAMS.h" +#include "SEAICE_SIZE.h" +#include "SEAICE_PARAMS.h" + +C !INPUT PARAMETERS: +C tracerId :: tracer identifier +C myThid :: my thread Id number + INTEGER tracerId + INTEGER myThid +CEOP + +C !LOCAL VARIABLES: + +C-- Set diagnostic suffix (4 character long) for the "tracerId" tracer + IF ( tracerId.EQ.GAD_HEFF ) THEN + SEAICE_DIAG_SUFX = 'HEFF' + ELSEIF( tracerId.EQ.GAD_AREA ) THEN + SEAICE_DIAG_SUFX = 'AREA' + ELSEIF( tracerId.EQ.GAD_SNOW ) THEN + SEAICE_DIAG_SUFX = 'SNOW' + ELSEIF( tracerId.EQ.GAD_SALT ) THEN + SEAICE_DIAG_SUFX = 'SSLT' + ELSE + SEAICE_DIAG_SUFX = 'aaaa' + ENDIF + + RETURN + END diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_model.F b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_model.F new file mode 100755 index 0000000..e516f07 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_model.F @@ -0,0 +1,390 @@ +#include "SEAICE_OPTIONS.h" +#ifdef ALLOW_AUTODIFF +# include "AUTODIFF_OPTIONS.h" +#endif + +CBOP +C !ROUTINE: SEAICE_MODEL + +C !INTERFACE: ========================================================== + SUBROUTINE SEAICE_MODEL( myTime, myIter, myThid ) + +C !DESCRIPTION: \bv +C *===========================================================* +C | SUBROUTINE SEAICE_MODEL | +C | o Time stepping of a dynamic/thermodynamic sea ice model. | +C | Dynamics solver: Zhang/Hibler, JGR, 102, 8691-8702, 1997 | +C | Thermodynamics: Hibler, MWR, 108, 1943-1973, 1980 | +C | Rheology: Hibler, JPO, 9, 815- 846, 1979 | +C | Snow: Zhang et al. , JPO, 28, 191- 217, 1998 | +C | Parallel forward ice model written by Jinlun Zhang PSC/UW| +C | & coupled into MITgcm by Dimitris Menemenlis (JPL) 2/2001| +C | zhang@apl.washington.edu / menemenlis@jpl.nasa.gov | +C *===========================================================* +C *===========================================================* + IMPLICIT NONE +C \ev + +C !USES: =============================================================== +#include "SIZE.h" +#include "EEPARAMS.h" +#include "DYNVARS.h" +#include "PARAMS.h" +#include "GRID.h" +#include "FFIELDS.h" +#include "SEAICE_SIZE.h" +#include "SEAICE_PARAMS.h" +#include "SEAICE.h" +#include "SEAICE_TRACER.h" +#ifdef ALLOW_EXF +# include "EXF_OPTIONS.h" +# include "EXF_FIELDS.h" +#endif +#ifdef ALLOW_AUTODIFF_TAMC +# include "tamc.h" +#endif +#ifdef HACK_FOR_GMAO_CPL +# include "SEAICE_LAYERS.h" +#endif /* HACK_FOR_GMAO_CPL */ + +C !INPUT PARAMETERS: =================================================== +C myTime :: Current time in simulation +C myIter :: Current iteration number in simulation +C myThid :: my Thread Id number + _RL myTime + INTEGER myIter + INTEGER myThid + +C !LOCAL VARIABLES: ==================================================== +C i,j,bi,bj :: Loop counters + INTEGER i, j + INTEGER bi, bj +#ifdef ALLOW_EXF + INTEGER grpDiag +#endif +#ifdef ALLOW_SITRACER + INTEGER iTr +#endif +#ifndef SEAICE_CGRID + _RL uLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL vLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) +#endif +CEOP + +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_ENTER( 'SEAICE_MODEL', myThid ) +#endif + +#ifdef ALLOW_THSICE + IF ( useThSice ) THEN +C-- Map thSice-variables to HEFF and AREA + CALL SEAICE_MAP_THSICE( myTime, myIter, myThid ) + ENDIF +#endif /* ALLOW_THSICE */ + +#ifdef ALLOW_EXF + IF ( useEXF ) THEN +C-- Winds are from pkg/exf, which does not update edges. + CALL EXCH_UV_AGRID_3D_RL( uwind, vwind, .TRUE., 1, myThid ) + IF ( useDiagnostics ) THEN +C- Fill-in EXF wind-stess diags, weighted by open-ocean fraction + grpDiag = -1 + IF ( SEAICEuseDYNAMICS ) grpDiag = 1 + CALL EXF_WEIGHT_SFX_DIAGS( + I AREA, grpDiag, myTime, myIter, myThid ) + ENDIF + ENDIF +#endif /* ALLOW_EXF */ + +#ifdef HACK_FOR_GMAO_CPL + CALL SEAICE_SAVE4GMAO( myTime, 1, myIter, myThid ) +#endif /* HACK_FOR_GMAO_CPL */ + +#ifdef ALLOW_AUTODIFF_TAMC + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + uIceNm1(i,j,bi,bj) = 0. _d 0 + vIceNm1(i,j,bi,bj) = 0. _d 0 +# ifdef ALLOW_SITRACER + DO iTr = 1, SItrMaxNum + SItrBucket(i,j,bi,bj,iTr) = 0. _d 0 + ENDDO +# endif + ENDDO + ENDDO + ENDDO + ENDDO +CADJ STORE uwind = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE vwind = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE heff = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte +# ifdef SEAICE_ALLOW_DYNAMICS +# ifdef SEAICE_CGRID +CADJ STORE hsnow = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE seaicemasku = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE seaicemaskv = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE fu = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE fv = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE uice = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE vice = comlev1, key=ikey_dynamics, kind=isbyte +cphCADJ STORE eta = comlev1, key=ikey_dynamics, kind=isbyte +cphCADJ STORE zeta = comlev1, key=ikey_dynamics, kind=isbyte +cph( +CADJ STORE dwatn = comlev1, key=ikey_dynamics, kind=isbyte +#ifdef SEAICE_ALLOW_BOTTOMDRAG +CADJ STORE cbotc = comlev1, key=ikey_dynamics, kind=isbyte +#endif /* SEAICE_ALLOW_BOTTOMDRAG */ +cccCADJ STORE press0 = comlev1, key=ikey_dynamics, kind=isbyte +cccCADJ STORE zmax = comlev1, key=ikey_dynamics, kind=isbyte +cccCADJ STORE zmin = comlev1, key=ikey_dynamics, kind=isbyte +cph) +# ifdef SEAICE_ALLOW_EVP +CADJ STORE seaice_sigma1 = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE seaice_sigma2 = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE seaice_sigma12 = comlev1, key=ikey_dynamics, kind=isbyte +# endif +# endif +# endif +# ifdef ALLOW_SITRACER +CADJ STORE siceload = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE sitracer = comlev1, key=ikey_dynamics, kind=isbyte +# endif +#endif /* ALLOW_AUTODIFF_TAMC */ + +C solve ice momentum equations and calculate ocean surface stress +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL( 'SEAICE_DYNSOLVER', myThid ) +#endif +#ifdef SEAICE_CGRID + CALL TIMER_START('SEAICE_DYNSOLVER [SEAICE_MODEL]',myThid) + CALL SEAICE_DYNSOLVER ( myTime, myIter, myThid ) + CALL TIMER_STOP ('SEAICE_DYNSOLVER [SEAICE_MODEL]',myThid) +#else + CALL TIMER_START('DYNSOLVER [SEAICE_MODEL]',myThid) + CALL DYNSOLVER ( myTime, myIter, myThid ) + CALL TIMER_STOP ('DYNSOLVER [SEAICE_MODEL]',myThid) +#endif /* SEAICE_CGRID */ + +C-- Apply ice velocity open boundary conditions +#ifdef ALLOW_OBCS +# ifndef DISABLE_SEAICE_OBCS + IF ( useOBCS ) CALL OBCS_ADJUST_UVICE( uice, vice, myThid ) +# endif /* DISABLE_SEAICE_OBCS */ +#endif /* ALLOW_OBCS */ + +#ifdef ALLOW_THSICE + IF ( useThSice ) THEN +#ifndef OLD_THSICE_CALL_SEQUENCE +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL( 'THSICE_DO_ADVECT', myThid ) +#endif + CALL THSICE_DO_ADVECT( 0, 0, myTime, myIter, myThid ) +#endif /* OLD_THSICE_CALL_SEQUENCE */ + ELSE +#endif /* ALLOW_THSICE */ +C-- Only call advection of heff, area, snow, and salt and +C-- growth for the generic 0-layer thermodynamics of seaice +C-- if useThSice=.false., otherwise the 3-layer Winton thermodynamics +C-- (called from DO_OCEANIC_PHYSICS) take care of this + +C NOW DO ADVECTION and DIFFUSION + IF ( SEAICEadvHeff .OR. SEAICEadvArea .OR. SEAICEadvSnow + & .OR. SEAICEadvSalt ) THEN +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL( 'SEAICE_ADVDIFF', myThid ) +#endif +#ifdef SEAICE_CGRID + CALL SEAICE_ADVDIFF( uIce, vIce, myTime, myIter, myThid ) +#else /* SEAICE_CGRID */ + CALL SEAICE_ADVDIFF( uLoc, vLoc, myTime, myIter, myThid ) +#endif /* SEAICE_CGRID */ + ENDIF + +C After advection, the sea ice variables may have unphysical values +C e.g., < 0, that are regularized here. Concentration as a special case +C may be > 1 in convergent motion and a ridging algorithm redistributes +C the ice to limit the concentration to 1. +#ifndef HACK_FOR_GMAO_CPL + CALL SEAICE_REG_RIDGE( myTime, myIter, myThid ) +#endif /* ndef HACK_FOR_GMAO_CPL */ + +#ifdef ALLOW_EXF + IF ( useEXF .AND. useDiagnostics ) THEN +C- Fill-in EXF surface flux diags, weighted by open-ocean fraction + grpDiag = -2 + IF ( usePW79thermodynamics ) grpDiag = 2 + CALL EXF_WEIGHT_SFX_DIAGS( + I AREA, grpDiag, myTime, myIter, myThid ) + ENDIF +#endif /* ALLOW_EXF */ + +#ifdef DISABLE_SEAICE_GROWTH + IF ( .TRUE. ) THEN +#else /* DISABLE_SEAICE_GROWTH */ +C thermodynamics growth +C must call growth after calling advection +C because of ugly time level business + IF ( usePW79thermodynamics ) THEN +#ifdef ALLOW_SEAICE_GROWTH_ADX +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL( 'SEAICE_GROWTH_ADX', myThid ) +#endif + CALL SEAICE_GROWTH_ADX( myTime, myIter, myThid ) +#else +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL( 'SEAICE_GROWTH', myThid ) +#endif + CALL SEAICE_GROWTH( myTime, myIter, myThid ) +#endif + ELSE +#endif /* DISABLE_SEAICE_GROWTH */ + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + DO j=1,sNy + DO i=1,sNx + sIceLoad(i,j,bi,bj) = HEFF(i,j,bi,bj)*SEAICE_rhoIce + & + HSNOW(i,j,bi,bj)*SEAICE_rhoSnow + ENDDO + ENDDO +c#ifdef SEAICE_CAP_ICELOAD +c sIceTooHeavy = rhoConst*drF(1) / 5. _d 0 +c DO j=1,sNy +c DO i=1,sNx +c sIceLoad(i,j,bi,bj) = MIN( sIceLoad(i,j,bi,bj), +c & sIceTooHeavy ) +c ENDDO +c ENDDO +c#endif + ENDDO + ENDDO + ENDIF + +#ifdef ALLOW_SITRACER +# ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE sitracer = comlev1, key=ikey_dynamics, kind=isbyte +# endif + CALL SEAICE_TRACER_PHYS ( myTime, myIter, myThid ) +#endif + +C-- Apply ice tracer open boundary conditions +#ifdef ALLOW_OBCS +# ifndef DISABLE_SEAICE_OBCS + IF ( useOBCS ) CALL OBCS_APPLY_SEAICE( myThid ) +# endif /* DISABLE_SEAICE_OBCS */ +#endif /* ALLOW_OBCS */ + +C-- Update overlap regions for a bunch of stuff + _EXCH_XY_RL( HEFF, myThid ) + _EXCH_XY_RL( AREA, myThid ) + _EXCH_XY_RL( HSNOW, myThid ) +#ifdef SEAICE_ITD + CALL EXCH_3D_RL( HEFFITD, nITD, myThid ) + CALL EXCH_3D_RL( AREAITD, nITD, myThid ) + CALL EXCH_3D_RL( HSNOWITD, nITD, myThid ) +#endif +#ifdef SEAICE_VARIABLE_SALINITY + _EXCH_XY_RL( HSALT, myThid ) +#endif +#ifdef ALLOW_SITRACER + DO iTr = 1, SItrNumInUse + _EXCH_XY_RL( SItracer(1-OLx,1-OLy,1,1,iTr),myThid ) + ENDDO +#endif + _EXCH_XY_RS(EmPmR, myThid ) + _EXCH_XY_RS(saltFlux, myThid ) + _EXCH_XY_RS(Qnet , myThid ) +#ifdef SHORTWAVE_HEATING + _EXCH_XY_RS(Qsw , myThid ) +#endif /* SHORTWAVE_HEATING */ +#ifdef ATMOSPHERIC_LOADING + IF ( useRealFreshWaterFlux ) + & _EXCH_XY_RS( sIceLoad, myThid ) +#endif + +#ifdef ALLOW_OBCS +C-- In case we use scheme with a large stencil that extends into overlap: +C no longer needed with the right masking in advection & diffusion S/R. +c IF ( useOBCS ) THEN +c DO bj=myByLo(myThid),myByHi(myThid) +c DO bi=myBxLo(myThid),myBxHi(myThid) +c CALL OBCS_COPY_TRACER( HEFF(1-OLx,1-OLy,bi,bj), +c I 1, bi, bj, myThid ) +c CALL OBCS_COPY_TRACER( AREA(1-OLx,1-OLy,bi,bj), +c I 1, bi, bj, myThid ) +c CALL OBCS_COPY_TRACER( HSNOW(1-OLx,1-OLy,bi,bj), +c I 1, bi, bj, myThid ) +#ifdef SEAICE_VARIABLE_SALINITY +c CALL OBCS_COPY_TRACER( HSALT(1-OLx,1-OLy,bi,bj), +c I 1, bi, bj, myThid ) +#endif +c ENDDO +c ENDDO +c ENDIF +#endif /* ALLOW_OBCS */ + +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics ) THEN +C diagnostics for "non-state variables" that are modified by +C the seaice model +# ifdef ALLOW_EXF +c CALL DIAGNOSTICS_FILL(UWIND ,'SIuwind ',0,1 ,0,1,1,myThid) +c CALL DIAGNOSTICS_FILL(VWIND ,'SIvwind ',0,1 ,0,1,1,myThid) +# endif +c CALL DIAGNOSTICS_FILL_RS(FU ,'SIfu ',0,1 ,0,1,1,myThid) +c CALL DIAGNOSTICS_FILL_RS(FV ,'SIfv ',0,1 ,0,1,1,myThid) + CALL DIAGNOSTICS_FILL_RS(EmPmR,'SIempmr ',0,1 ,0,1,1,myThid) + CALL DIAGNOSTICS_FILL_RS(Qnet ,'SIqnet ',0,1 ,0,1,1,myThid) + CALL DIAGNOSTICS_FILL_RS(Qsw ,'SIqsw ',0,1 ,0,1,1,myThid) + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ + +#ifdef ALLOW_THSICE +C endif .not.useThSice + ENDIF +#endif /* ALLOW_THSICE */ +CML This has already been done in seaice_ocean_stress/ostres, so why repeat? +CML CALL EXCH_UV_XY_RS(fu,fv,.TRUE.,myThid) + +#ifdef HACK_FOR_GMAO_CPL +C- Save advective increments: + CALL SEAICE_SAVE4GMAO( myTime, 2, myIter, myThid ) +C- for now, reset increments to zero +c CALL SEAICE_SAVE4GMAO( myTime, 0, myIter, myThid ) + +C- Compute SI_FRZMLT, the available heat (W/m^2) to freeze (>0) or melt (<0) +C sea ice so that surface level ocean reaches freezing temperature. +C FRZMLT = (-0.054*SW - TW) * (MAPL_RHO_SEAWATER*MAPL_CAPWTR*10.)/DT +C where -0.054*SW is the CICE freezing point of sea water in deg C. +C The freezing potential (>0) is subtracted from QNET. + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + SI_FRZMLT(i,j,bi,bj) = maskC(i,j,1,bi,bj) * + & (-0.054*salt(i,j,1,bi,bj) - theta(i,j,1,bi,bj)) * + & rhoNil * HeatCapacity_Cp * drF(1) / deltaT + Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - + & MAX( SI_FRZMLT(i,j,bi,bj) , zeroRL ) + ENDDO + ENDDO + ENDDO + ENDDO + +#endif /* HACK_FOR_GMAO_CPL */ + +#ifdef ALLOW_EXF +# ifdef ALLOW_AUTODIFF_TAMC +# if (defined (ALLOW_AUTODIFF_MONITOR)) + CALL EXF_ADJOINT_SNAPSHOTS( 3, myTime, myIter, myThid ) +# endif +# endif +#endif + +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_LEAVE( 'SEAICE_MODEL', myThid ) +#endif + + RETURN + END diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_save4gmao.F b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_save4gmao.F new file mode 100755 index 0000000..bd61105 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/code/seaice_save4gmao.F @@ -0,0 +1,261 @@ +#include "SEAICE_OPTIONS.h" + +CBOP +C !ROUTINE: SEAICE_SAVE4GMAO +C !INTERFACE: + SUBROUTINE SEAICE_SAVE4GMAO( + I myTime, seqFlag, myIter, myThid ) + +C !DESCRIPTION: \bv +C *==========================================================* +C | S/R SEAICE_SAVE4GMAO +C | o Save Seaice Advective Increment for GMAO Coupling +C *==========================================================* +C \ev + +C !USES: + IMPLICIT NONE + +C == Global variables === +#include "SIZE.h" +#include "EEPARAMS.h" +#include "PARAMS.h" +#include "GRID.h" +#ifdef ALLOW_SEAICE +# include "SEAICE_SIZE.h" +# include "SEAICE_PARAMS.h" +# include "SEAICE.h" +# ifdef HACK_FOR_GMAO_CPL +# include "SEAICE_LAYERS.h" +# endif +#endif /* ALLOW_SEAICE */ + +C !INPUT/OUTPUT PARAMETERS: +C myTime :: Current time of simulation ( s ) +C seqFlag :: flag that indicate where this S/R is called from: +C :: =0 reset inc. to zero (from initialise_varia) +C :: =1 called from the beginning of SEAICE_MODEL +C :: =2 called from the end of SEAICE_MODEL +C myIter :: Iteration number +C myThid :: my Thread Id number + _RL myTime + INTEGER seqFlag + INTEGER myIter + INTEGER myThid +CEOP + +#ifdef HACK_FOR_GMAO_CPL +# ifdef ALLOW_SEAICE +C !LOCAL VARIABLES: +C i,j,bi,bj :: Loop counters + INTEGER i, j, bi, bj + INTEGER l, n +c _RL tmpFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy) +c _RL fld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| + + IF ( seqFlag.EQ.0 ) THEN +C-- Initialise to zero Advective Increments + + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + + DO n=1,nITD +C- start loop on ice-category "n" + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + SIadv_Area (i,j,n,bi,bj) = 0. _d 0 + SIadv_Heff (i,j,n,bi,bj) = 0. _d 0 + SIadv_Hsnow (i,j,n,bi,bj) = 0. _d 0 + SIadv_meltPd(i,j,n,bi,bj) = 0. _d 0 + SIadv_iceAge(i,j,n,bi,bj) = 0. _d 0 + SIadv_tIces (i,j,n,bi,bj) = 0. _d 0 + ENDDO + ENDDO + DO l=1,nIceLayers + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + SIadv_qIce (i,j,l,n,bi,bj) = 0. _d 0 + ENDDO + ENDDO + ENDDO + DO l=1,nSnowLayers + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + SIadv_qSnow (i,j,l,n,bi,bj) = 0. _d 0 + ENDDO + ENDDO + ENDDO +C- end loop on ice-category "n" + ENDDO + + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + SIadv_skinS (i,j,bi,bj) = 0. _d 0 + SIadv_skinH (i,j,bi,bj) = 0. _d 0 + ENDDO + ENDDO + +C- end bi,bj loops + ENDDO + ENDDO + + ENDIF + +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| + + IF ( seqFlag.EQ.1 ) THEN +C-- Save seaice fields as entering seaice model + + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + + DO n=1,nITD +C- start loop on ice-category "n" + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + SIadv_Area (i,j,n,bi,bj) = AREAITD (i,j,n,bi,bj) + SIadv_Heff (i,j,n,bi,bj) = HEFFITD (i,j,n,bi,bj) + SIadv_Hsnow (i,j,n,bi,bj) = HSNOWITD(i,j,n,bi,bj) + SIadv_meltPd(i,j,n,bi,bj) = SImeltPd(i,j,n,bi,bj) + SIadv_iceAge(i,j,n,bi,bj) = SIiceAge(i,j,n,bi,bj) + SIadv_tIces (i,j,n,bi,bj) = TICES (i,j,n,bi,bj) + ENDDO + ENDDO + DO l=1,nIceLayers + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + SIadv_qIce (i,j,l,n,bi,bj) = SIqIce (i,j,l,n,bi,bj) + ENDDO + ENDDO + ENDDO + DO l=1,nSnowLayers + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + SIadv_qSnow(i,j,l,n,bi,bj) = SIqSnow(i,j,l,n,bi,bj) + ENDDO + ENDDO + ENDDO +C- end loop on ice-category "n" + ENDDO + +c DO j=1-OLy,sNy+OLy +c DO i=1-OLx,sNx+OLx +c SIadv_skinS (i,j,bi,bj) = SIskinS(i,j,bi,bj) +c SIadv_skinH (i,j,bi,bj) = SIskinH(i,j,bi,bj) +c ENDDO +c ENDDO + +C- end bi,bj loops + ENDDO + ENDDO + +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics ) THEN + CALL DIAGNOSTICS_FILL( oceWeight,'CPLoWGHT',0,1,0,1,1,myThid ) + CALL DIAGNOSTICS_FILL( TICES, 'SItIces ', + I 0, nITD, 0,1,1, myThid ) + n = nIceLayers*nITD + CALL DIAGNOSTICS_FILL( SIqIce, 'SIqIce ', + I 0, n , 0,1,1, myThid ) + n = nSnowLayers*nITD + CALL DIAGNOSTICS_FILL( SIqSnow, 'SIqSnow ', + I 0, n , 0,1,1, myThid ) + CALL DIAGNOSTICS_FILL( SImeltPd, 'SImeltPd', + I 0, nITD, 0,1,1, myThid ) + CALL DIAGNOSTICS_FILL( SIiceAge, 'SIiceAge', + I 0, nITD, 0,1,1, myThid ) + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ + + ENDIF + +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| + + IF ( seqFlag.EQ.2 ) THEN +C-- Calculate Advective Increments + + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + + DO n=1,nITD +C- start loop on ice-category "n" + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + SIadv_Area (i,j,n,bi,bj) = AREAITD (i,j,n,bi,bj) + & - SIadv_Area (i,j,n,bi,bj) + SIadv_Heff (i,j,n,bi,bj) = HEFFITD (i,j,n,bi,bj) + & - SIadv_Heff (i,j,n,bi,bj) + SIadv_Hsnow (i,j,n,bi,bj) = HSNOWITD(i,j,n,bi,bj) + & - SIadv_Hsnow (i,j,n,bi,bj) + SIadv_meltPd(i,j,n,bi,bj) = SImeltPd(i,j,n,bi,bj) + & - SIadv_meltPd(i,j,n,bi,bj) + SIadv_iceAge(i,j,n,bi,bj) = SIiceAge(i,j,n,bi,bj) + & - SIadv_iceAge(i,j,n,bi,bj) + SIadv_tIces (i,j,n,bi,bj) = TICES (i,j,n,bi,bj) + & - SIadv_tIces (i,j,n,bi,bj) + ENDDO + ENDDO + DO l=1,nIceLayers + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + SIadv_qIce (i,j,l,n,bi,bj) = SIqIce (i,j,l,n,bi,bj) + & - SIadv_qIce (i,j,l,n,bi,bj) + ENDDO + ENDDO + ENDDO + DO l=1,nSnowLayers + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + SIadv_qSnow(i,j,l,n,bi,bj) = SIqSnow(i,j,l,n,bi,bj) + & - SIadv_qSnow(i,j,l,n,bi,bj) + ENDDO + ENDDO + ENDDO +C- end loop on ice-category "n" + ENDDO + +c DO j=1-OLy,sNy+OLy +c DO i=1-OLx,sNx+OLx +c SIadv_skinS(i,j,bi,bj) = SIskinS(i,j,bi,bj) +c & - SIadv_skinS(i,j,bi,bj) +c SIadv_skinH(i,j,bi,bj) = SIskinH(i,j,bi,bj) +c & - SIadv_skinH(i,j,bi,bj) +c ENDDO +c ENDDO + +C- end bi,bj loops + ENDDO + ENDDO + +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics ) THEN + CALL DIAGNOSTICS_FILL( SIadv_Area, 'SI_dArea', + I 0, nITD, 0,1,1, myThid ) + CALL DIAGNOSTICS_FILL( SIadv_Heff, 'SI_dHeff', + I 0, nITD, 0,1,1, myThid ) + CALL DIAGNOSTICS_FILL( SIadv_Hsnow, 'SI_dHsnw', + I 0, nITD, 0,1,1, myThid ) + CALL DIAGNOSTICS_FILL( SIadv_tIces, 'SI_dTIce', + I 0, nITD, 0,1,1, myThid ) + n = nIceLayers*nITD + CALL DIAGNOSTICS_FILL( SIadv_qIce, 'SI_dQIce', + I 0, n , 0,1,1, myThid ) + n = nSnowLayers*nITD + CALL DIAGNOSTICS_FILL( SIadv_qSnow, 'SI_dQSnw', + I 0, n , 0,1,1, myThid ) + CALL DIAGNOSTICS_FILL( SIadv_meltPd, 'SI_dMPnd', + I 0, nITD, 0,1,1, myThid ) + CALL DIAGNOSTICS_FILL( SIadv_iceAge, 'SI_dIcAg', + I 0, nITD, 0,1,1, myThid ) + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ + + ENDIF + +# endif /* ALLOW_SEAICE */ +#endif /* HACK_FOR_GMAO_CPL */ + + RETURN + END diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data new file mode 100755 index 0000000..31abce9 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data @@ -0,0 +1,132 @@ +# ==================== +# | Model parameters | +# ==================== +# +# Continuous equation parameters + &PARM01 + tRef = 3*23.,3*22.,21.,2*20.,19.,2*18.,17.,2*16.,15.,14.,13., + 12.,11.,2*9.,8.,7.,2*6.,2*5.,3*4.,3*3.,4*2.,12*1., + sRef = 50*34.5, + no_slip_sides = .TRUE., + no_slip_bottom = .TRUE., +# + viscAr=0.5E-4, +# + viscAh=1.E0, + viscAhGrid=1.E-2, +# viscAh=2.0e4, +# + diffKhT=1.E1, + diffKrT=1.E-5, + diffKhS=1.E1, + diffKrS=1.E-5, +# +### diffKrBL79surf=0.1E-4, +### diffKrBL79deep=1.0E-4, + bottomDragQuadratic = 0.001, +#when using ggl90 + ivdc_kappa=10., + implicitDiffusion=.TRUE., + implicitViscosity=.TRUE., + viscC4Leith=1.5, + viscC4Leithd=1.5, + viscA4GridMax=0.5, + useAreaViscLength=.TRUE., + useRealFreshWaterFlux=.TRUE., +# balanceThetaClimRelax=.TRUE., +# balanceSaltClimRelax=.TRUE., +# balanceEmPmR=.TRUE., +# balanceQnet=.TRUE., + allowFreezing=.FALSE., +### hFacInf=0.2, +### hFacSup=2.0, +# hFacMin=.2, +# hFacMinDr=5., + hFacMinDr=50., + hFacMin=0.3, + hFacInf=0.1, + hFacSup=5., + select_rStar=2, + nonlinFreeSurf=4, + gravity=9.81, + rhonil=1029., + rhoConst=1029., + rhoConstFresh=1000., + convertFW2Salt=-1., + eosType='JMD95Z', + implicitFreeSurface=.TRUE., + exactConserv=.TRUE., + useSingleCpuIO=.TRUE., + tempAdvScheme=30, + saltAdvScheme=30, + tempVertAdvScheme=3, + saltVertAdvScheme=3, + tempImplVertAdv=.TRUE., + saltImplVertAdv=.TRUE., + staggerTimeStep=.TRUE., + vectorInvariantMomentum=.TRUE., +#when using the cd scheme: +# useCDscheme=.TRUE., + useJamartWetPoints=.TRUE., + readBinaryPrec=32, + writeBinaryPrec=32, + debugLevel=1, + / + +# Elliptic solver parameters + &PARM02 + cg2dMaxIters=300, +#cg2dTargetResWunit=1.E-12, + / + +# Time stepping parameters + &PARM03 + nIter0=2, +#26yr + nTimeSteps=227902, + forcing_In_AB=.FALSE., + momDissip_In_AB=.FALSE., +#when using the cd scheme: +# epsAB_CD = 0.25, +# tauCD=172800.0, + deltaT = 1200., +#when using ab2: + abEps = 0.1, +#when using ab3: +# doAB_onGtGs=.FALSE., +# alph_AB=0.5, +# beta_AB=0.281105, +# + pChkptFreq =31536000.0, + chkptFreq =31536000.0, +# taveFreq =31536000.0, +# dumpFreq =31536000.0, + monitorFreq=3600, + dumpInitAndLast = .TRUE., + pickupStrictlyMatch=.FALSE., + / + +# Gridding parameters + &PARM04 + usingCurvilinearGrid=.TRUE., + delR = + 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, + 10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04, 19.82, 24.85, + 31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, + 93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, + 139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, + 341.50,364.50,387.50,410.50,433.50,456.50, + / + +# Input datasets + &PARM05 + diffKrFile='ecco_v5mg_diffkr.data', +#adTapeDir='tapes', +#bathyFile ='bathy_eccollc_90x50.bin', + bathyFile ='bathy270_filled_noCaspian_r4', + hydrogThetaFile='T_OWPv1_M_eccollc_90x50.bin', + hydrogSaltFile ='S_OWPv1_M_eccollc_90x50.bin', +# viscA4Dfile ='viscA4Dfld_eccollc_90x50.bin', +# viscA4Zfile ='viscA4Zfld_eccollc_90x50.bin', +# + / diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.diagnostics b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.diagnostics new file mode 100755 index 0000000..07c5045 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.diagnostics @@ -0,0 +1,127 @@ +# Diagnostic Package Choices +#-------------------- +# dumpAtLast (logical): always write output at the end of simulation (default=F) +# diag_mnc (logical): write to NetCDF files (default=useMNC) +#--for each output-stream: +# fileName(n) : prefix of the output file name (max 80c long) for outp.stream n +# frequency(n):< 0 : write snap-shot output every |frequency| seconds +# > 0 : write time-average output every frequency seconds +# timePhase(n) : write at time = timePhase + multiple of |frequency| +# averagingFreq : frequency (in s) for periodic averaging interval +# averagingPhase : phase (in s) for periodic averaging interval +# repeatCycle : number of averaging intervals in 1 cycle +# levels(:,n) : list of levels to write to file (Notes: declared as REAL) +# when this entry is missing, select all common levels of this list +# fields(:,n) : list of selected diagnostics fields (8.c) in outp.stream n +# (see "available_diagnostics.log" file for the full list of diags) +# missing_value(n) : missing value for real-type fields in output file "n" +# fileFlags(n) : specific code (8c string) for output file "n" +#-------------------- +# oceTAUX zonal surface wind stress, >0 increases uVel (N/m^2) +# oceTAUY meridional surf. wind stress, >0 increases vVel (N/m^2) +# oceFWflx net surface Fresh-Water flux into ocean, >0 decreases salinity (kg/m^2/s) +# oceSflux net surface Salt flux into the ocean, >0 increases salinity (g/m^2/s) +# oceQnet net surface heat flux into the ocean, >0 increases theta (W/m^2) +# oceQsw net Short-Wave radiation (+=down), >0 increases theta (W/m^2) +# SSS Sea Surface Salinity (g/kg) +# SST Sea Surface Temperature (degC) +# UVEL1 Zonal Surface Velocity (m/s) +# VVEL1 Meridional Surface Velocity (m/s) + + &DIAGNOSTICS_LIST + + frequency(1) = 86400., + fields(1:8,1) = 'oceTAUX ','oceTAUY ','oceFWflx','oceSflux','oceQnet ','oceQsw ','ETAN ','atmPload', + filename(1) = 'state_2d_set1', + timePhase(1) = 27000., + + frequency(2) = 86400., + fields(1:7,2) = 'SALT ','THETA ','UVELMASS','VVELMASS','WVELMASS','GM_PsiX ','GM_PsiY ', + filename(2) = 'state_3d_set1', + timePhase(2) = 27000., + + frequency(3) = -86400., + fields(1:2,3) = 'SALT ','THETA ', + filename(3) = 'snap_3d_set1', + timePhase(3) = 27000., + + frequency(4) = -86400.0, + fields(1,4) = 'ETAN ', + filename(4) = 'snap_2d_set2', + timePhase(4) = 27000., + + frequency(5) = 86400.0, + fields(1:8,5) = 'DFxE_TH ','DFyE_TH ','ADVx_TH ','ADVy_TH ', + 'DFxE_SLT','DFyE_SLT','ADVx_SLT','ADVy_SLT', + filename(5) = 'state_3d_set2', + timePhase(5) = 27000., + + frequency(6) = 86400.0, + fields(1:7,6) = 'ADVr_TH ','DFrE_TH ','DFrI_TH ', + 'ADVr_SLT','DFrE_SLT','DFrI_SLT', + filename(6) = 'state_3d_set4', + timePhase(6) = 27000., + + frequency(7) = 86400.0, + fields(1:14,7) = 'SIarea ','SIheff ','SIhsnow ','sIceLoad', + 'SIuice ','SIvice ','SItaux ','SItauy ', + 'ADVxHEFF','ADVyHEFF', + 'ADVxAREA','ADVyAREA', + 'ADVxSNOW','ADVySNOW', + fileName(7) = 'iceDiag', + missing_value(7) = -999., + timePhase(7) = 27000., + + frequency(8) = -86400., + fields(1:9,8) = 'SIarea ','SIheff ','SIhsnow ','sIceLoad', + 'SIuice ','SIvice ','SItaux ','SItauy ', + 'CPLoWGHT', + missing_value(8) = -999., +# fileName(8) = 'iceInst', + timePhase(8) = 27000., + + frequency(9) = 86400., + fields(1:11,9) = 'SIareaN ','SIheffN ','SIhsnowN','SIqSnow ', + 'SImeltPd','SItIces ','SIiceAge', + 'SI_dArea','SI_dHeff','SI_dHsnw','SI_dQSnw', +# 'SI_dMPnd', +# 'SI_dTIce','SI_dIcAg', + timePhase(9) = 0., + fileFlags(9) = 'D ', +# fileName(9) = 'iceNcat', + timePhase(9) = 27000., + + frequency(10) = 86400., + fields(1:2,10) = 'SIqIce ', + 'SI_dQIce', + timePhase(10) = 0., +# fileName(10) = 'iceEnerg', + timePhase(10) = 27000., + + frequency(11) = 86400., + fields(1:2,11) = 'SIqSnow ', + 'SI_dQSnw', + timePhase(11) = 0., +# fileName(11) = 'snowEnerg', + timePhase(11) = 27000., + + & + + &DIAG_STATIS_PARMS + stat_fields(1:6,1) = 'ETAN ', 'THETA ','SALT ','oceQnet ','oceFWflx','oceSflux', + stat_fName(1) = 'dynStDiag', + stat_freq(1)= 86400.0, + stat_phase(2) = 27000., + + stat_fields(1:21,2) = 'SIarea ','SIheff ','SIhsnow ','sIceLoad', + 'SIuice ','SIvice ','SItaux ','SItauy ', + 'SIareaN ','SIheffN ','SIhsnowN', + 'SIqSnow ','SIqIce ','SImeltPd', + 'SItIces ','SIiceAge', + 'SI_dArea','SI_dHeff','SI_dHsnw', + 'SI_dQSnw','SI_dQIce', +# 'SI_dMPnd','SI_dTIce','SI_dIcAg', + stat_fName(2) = 'iceStDiag', + stat_freq(2) = 86400., + stat_phase(2) = 27000., + & diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.exch2 b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.exch2 new file mode 100755 index 0000000..5fb4d72 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.exch2 @@ -0,0 +1,318 @@ +# EXCH2 Package: Wrapper-2 User Choice +#-------------------- +# preDefTopol :: pre-defined Topology selector: +# :: = 0 : topology defined from processing "data.exch2"; +# :: = 1 : simple, single facet topology; +# :: = 2 : customized topology (w2_set_myown_facets) +# :: = 3 : 6-facet Cube (3 face-dims: nRed, nGreen, nBlue). +# dimsFacets :: facet pair of dimensions (n1x,n1y, n2x,n2y ...) +# facetEdgeLink :: Face-Edge connectivity map: +# facetEdgeLink(i,j)=XX.1 : face(j)-edge(i) (i=1,2,3,4 <==> N,S,E,W) +# is connected to Northern edge of face "XX" ; similarly, +# = XX.2 : to Southern.E, XX.3 = Eastern.E, XX.4 = Western.E of face "XX" +# blankList :: List of "blank" tiles +# W2_mapIO :: global map IO selector (-1 = old type ; 0 = 1 long line in X +# :: 1 = compact, mostly in Y dir) +# W2_printMsg :: option for information messages printing +# :: <0 : write to log file ; =0 : minimum print ; +# :: =1 : no duplicated print ; =2 : all processes do print +#-------------------- + &W2_EXCH2_PARM01 + W2_printMsg= 0, + W2_mapIO = 1, +# + preDefTopol=0, +#-- 5 facets llc_270 topology (drop facet 6 and its connection): + dimsFacets = 270, 810, 270, 810, 270, 270, 810, 270, 810, 270, 0, 0, + facetEdgeLink(1,1)= 3.4, 0. , 2.4, 5.1, + facetEdgeLink(1,2)= 3.2, 0. , 4.2, 1.3, + facetEdgeLink(1,3)= 5.4, 2.1, 4.4, 1.1, + facetEdgeLink(1,4)= 5.2, 2.3, 0. , 3.3, + facetEdgeLink(1,5)= 1.4, 4.1, 0. , 3.1, +# + blankList = 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 11, + 12, + 13, + 14, + 15, + 16, + 17, + 18, + 19, + 20, + 21, + 22, + 23, + 24, + 25, + 26, + 27, + 29, + 30, + 31, + 32, + 33, + 34, + 35, + 36, + 39, + 40, + 41, + 42, + 43, + 44, + 45, + 54, + 133, + 142, + 150, + 151, + 159, + 160, + 161, + 166, + 167, + 168, + 169, + 175, + 176, + 178, + 189, + 198, + 207, + 214, + 215, + 216, + 224, + 225, + 234, + 244, + 245, + 246, + 247, + 248, + 249, + 250, + 251, + 252, + 253, + 254, + 255, + 256, + 257, + 258, + 259, + 260, + 261, + 262, + 263, + 264, + 265, + 266, + 267, + 268, + 269, + 270, + 271, + 272, + 273, + 274, + 275, + 276, + 277, + 278, + 279, + 280, + 281, + 282, + 283, + 284, + 285, + 286, + 287, + 288, + 289, + 292, + 293, + 294, + 295, + 296, + 297, + 368, + 369, + 417, + 418, + 419, + 420, + 424, + 425, + 426, + 427, + 428, + 429, + 433, + 434, + 435, + 436, + 437, + 438, + 442, + 443, + 444, + 445, + 446, + 447, + 448, + 449, + 451, + 452, + 453, + 454, + 455, + 456, + 457, + 458, + 460, + 461, + 462, + 463, + 464, + 465, + 466, + 467, + 469, + 470, + 471, + 472, + 473, + 474, + 475, + 476, + 477, + 478, + 481, + 482, + 483, + 484, + 485, + 486, + 495, + 543, + 551, + 559, + 568, + 569, + 589, + 590, + 591, + 592, + 593, + 594, + 595, + 596, + 617, + 618, + 619, + 620, + 621, + 645, + 646, + 647, + 648, + 673, + 674, + 675, + 700, + 701, + 702, + 727, + 728, + 729, + 731, + 754, + 755, + 756, + 758, + 780, + 781, + 782, + 783, + 784, + 785, + 807, + 808, + 809, + 810, + 811, + 812, + 813, + 814, + 834, + 835, + 836, + 837, + 839, + 840, + 841, + 842, + 843, + 861, + 862, + 863, + 864, + 866, + 867, + 868, + 869, + 870, + 871, + 888, + 889, + 890, + 891, + 896, + 897, + 898, + 914, + 915, + 916, + 917, + 918, + 924, + 942, + 943, + 944, + 945, + 968, + 969, + 970, + 971, + 972, + 983, + 984, + 985, + 986, + 987, + 996, + 997, + 998, + 999, + 1011, + 1012, + 1023, + 1024, + 1025, + 1026, + 1051, + 1052, + & diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.ggl90 b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.ggl90 new file mode 100755 index 0000000..59b518c --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.ggl90 @@ -0,0 +1,15 @@ +# ===================================================================== +# | Parameters for Gaspar et al. (1990)'s TKE vertical mixing scheme | +# ===================================================================== + &GGL90_PARM01 +# GGL90taveFreq = 345600000., +# GGL90dumpFreq = 86400., +# GGL90writeState=.FALSE., +# GGL90diffTKEh=3.e3, + GGL90alpha=30., +# GGL90TKEFile = 'TKE_init.bin', + GGL90TKEmin = 1.e-7, + GGL90TKEbottom = 1.e-6, + mxlMaxFlag =2, + mxlSurfFlag=.TRUE., + / diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.gmredi b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.gmredi new file mode 100755 index 0000000..8cd0303 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.gmredi @@ -0,0 +1,38 @@ +# GM+Redi package parameters: +# GM_Small_Number :: epsilon used in computing the slope +# GM_slopeSqCutoff :: slope^2 cut-off value + +#-from MOM : +# GM_background_K: G & Mc.W diffusion coefficient +# GM_maxSlope : max slope of isopycnals +# GM_Scrit : transition for scaling diffusion coefficient +# GM_Sd : half width scaling for diffusion coefficient +# GM_taper_scheme: slope clipping or one of the tapering schemes +# GM_Kmin_horiz : horizontal diffusion minimum value + +#-Option parameters (needs to "define" options in GMREDI_OPTIONS.h") +# GM_isopycK : isopycnal diffusion coefficient (default=GM_background_K) +# GM_AdvForm : turn on GM Advective form (default=Skew flux form) + + &GM_PARM01 + GM_Small_Number = 1.D-20, + GM_slopeSqCutoff = 1.D+08, + GM_AdvForm = .TRUE., + GM_isopycK = 100., + GM_background_K = 100., + GM_maxSlope = 4.D-3, + GM_taper_scheme = 'stableGmAdjTap', + GM_Kmin_horiz = 10., + GM_Scrit = 4.D-3, + GM_Sd = 1.D-3, + GM_K3dGMFile = 'ecco_v5mg_kapgm.data', + GM_K3dRediFile = 'ecco_v5mg_kapredi.data', +# +### GM_Visbeck_alpha = 1.5D-2, +### GM_Visbeck_alpha = 0.D0, +### GM_Visbeck_length = 2.D+5, +### GM_Visbeck_depth = 1.D+3, +### GM_Visbeck_maxval_K= 2.5D+3, + / + + diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.pkg b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.pkg new file mode 100755 index 0000000..2805dc1 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.pkg @@ -0,0 +1,8 @@ +# Packages + &PACKAGES + useGMRedi = .TRUE., + useGGL90 = .TRUE., + useSALT_PLUME = .TRUE., + useDiagnostics = .TRUE., + useSEAICE = .TRUE., + & diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.salt_plume b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.salt_plume new file mode 100755 index 0000000..8cb64c5 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.salt_plume @@ -0,0 +1,6 @@ + &SALT_PLUME_PARM01 +# SaltPlumeCriterion = 0.4D0, + SPsalFRAC= 0.5D0, +#SPsalFRAC= 0.25D0, +#SPsalFRAC= 0.0D0, + / diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.seaice b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.seaice new file mode 100755 index 0000000..7c8e416 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/data.seaice @@ -0,0 +1,50 @@ +# SEAICE parameters + &SEAICE_PARM01 + usePW79thermodynamics= .FALSE., + + SEAICEpresH0=2., + SEAICEpresPow0=1, + SEAICEpresPow1=1, + + SEAICE_strength = 2.25e4, + + SEAICE_no_slip = .TRUE., + + SEAICE_drag=0.001, + OCEAN_drag=0.001, + +##### + SEAICEuseTILT=.FALSE., + SEAICE_multDim=1, + SEAICErestoreUnderIce=.TRUE., + + LSR_ERROR = 2.e-4, + SEAICEuseDYNAMICS = .TRUE., + SEAICEadvScheme = 33, + SEAICEuseFluxForm = .TRUE., + SEAICEadvSnow = .TRUE., + SEAICEdiffKhHeff = 400., + SEAICEdiffKhArea = 400., + SEAICEdiffKhSnow = 400., + SEAICEuseFlooding = .TRUE., + SEAICE_mcPheePiston= 3.858024691358025E-05, + SEAICE_frazilFrac = 1., + SEAICE_mcPheeTaper = 0., + SEAICE_areaLossFormula=2, + SEAICEheatConsFix = .TRUE., + SEAICE_tempFrz0 = -1.96, + SEAICE_dTempFrz_dS = 0., + SEAICEuseMetricTerms = .TRUE., + +#changes needed due to pkg/seaice changes +#from checkpoint67c to checkpoint67d + SEAICEscaleSurfStress=.FALSE., + SEAICEaddSnowMass=.FALSE., + SEAICE_OLx=0, + SEAICE_OLy=0, + SEAICEetaZmethod=0, + SEAICE_waterDrag=5.34499514091350826D-3, + & + + &SEAICE_PARM03 + & diff --git a/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/eedata b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/eedata new file mode 100755 index 0000000..0c1e412 --- /dev/null +++ b/MIT_GEOS5PlugMod/configs/c180_llc270_01/input/eedata @@ -0,0 +1,13 @@ +# Example "eedata" file +# Lines beginning "#" are comments +# nTx - No. threads per process in X +# nTy - No. threads per process in Y + &EEPARMS + useCubedSphereExchange=.TRUE., + nTx=1, + nTy=1, +#debugMode=.TRUE., + / +# Note: Some systems use & as the +# namelist terminator. Other systems +# use a / character (as shown here).