From 45315406aaced376d8a901f197adaa94c702ec86 Mon Sep 17 00:00:00 2001 From: Martin Losch <30285667+mjlosch@users.noreply.github.com> Date: Thu, 3 Aug 2023 18:50:12 +0200 Subject: [PATCH] replace flag SEAICE_ALLOW_DYNAMICS by more consistent SEAICE_CGRID and SEAICE_BGRID_DYNAMICS (#734) * compile code only if SEAICE_ALLOW_DYNAMCS is defined - also refine usage of cpp flags * avoid unused variable warning * define and initialise variables only with SEAICE_ALLOW_DYNAMICS def. - avoid unnecessary memory overhead - re-organise common blocks and initialisation to reduce number cpp flags * add new flag SEAICE_BGRID_DYNAMICS, remove SEAICE_ALLOW_DYNAMICS - replace ifdef SEAICE_ALLOW_DYNAMICS with if ( define SEAICE_CGRID || defined SEAICE_BGRID_DYNAMICS ) - problem: cannot add check if SEAICE_CGRID and SEAICE_BGRID_DYNAMICS are both defined, because in this case, the code does not (yet) compile * amending previous incomplete commit * replace U/VICEC by u/vIceB for B-grid code avoid conflicts with C-grid code * more adjustments to make it work with neither B or C-grid defined and AD-ocde * fix a comment * update SEAICE_OPTIONS.h files in verification experiments to be as similar as possible to currrent pkg/seaice/SEAICE_OPTIONS.h * minor: some spacing to help read the code * avoid un-used variables * fix previous commit * minor cleaning * add comments and exclude more diagnostics code ifnde SEAICE_CGRID - explain potentially unclear CPP-flag logic - add B-grid/C-grid version comment and clean up header * fix a typo in comments * Document adding SEAICE_BGRID_DYNAMICS option --------- Co-authored-by: Jean-Michel Campin --- doc/tag-index | 6 + model/src/the_main_loop.F | 2 +- pkg/autodiff/autodiff_ini_model_io.F | 2 +- pkg/seaice/SEAICE.h | 142 ++++++------ pkg/seaice/SEAICE_OPTIONS.h | 27 ++- pkg/seaice/dynsolver.F | 21 +- pkg/seaice/lsr.F | 82 ++++--- pkg/seaice/ostres.F | 16 +- pkg/seaice/seaice_ad_check_lev1_dir.h | 10 +- pkg/seaice/seaice_ad_check_lev2_dir.h | 4 +- pkg/seaice/seaice_ad_check_lev3_dir.h | 4 +- pkg/seaice/seaice_ad_check_lev4_dir.h | 4 +- pkg/seaice/seaice_ad_dump.F | 8 +- pkg/seaice/seaice_advdiff.F | 30 +-- pkg/seaice/seaice_calc_ice_strength.F | 2 + pkg/seaice/seaice_calc_lhs.F | 5 +- pkg/seaice/seaice_calc_residual.F | 4 +- pkg/seaice/seaice_calc_rhs.F | 5 +- pkg/seaice/seaice_calc_strainrates.F | 2 - pkg/seaice/seaice_calc_viscosities.F | 4 +- pkg/seaice/seaice_check.F | 29 ++- pkg/seaice/seaice_dynsolver.F | 111 +++++---- pkg/seaice/seaice_evp.F | 6 +- pkg/seaice/seaice_freedrift.F | 6 +- pkg/seaice/seaice_get_dynforcing.F | 3 - pkg/seaice/seaice_init_fixed.F | 65 +++--- pkg/seaice/seaice_init_varia.F | 131 ++++++----- pkg/seaice/seaice_jfnk.F | 4 +- pkg/seaice/seaice_krylov.F | 4 +- pkg/seaice/seaice_lsr.F | 13 -- pkg/seaice/seaice_model.F | 54 ++--- pkg/seaice/seaice_monitor.F | 10 +- pkg/seaice/seaice_monitor_ad.F | 10 +- pkg/seaice/seaice_ocean_stress.F | 7 +- pkg/seaice/seaice_oceandrag_coeffs.F | 4 +- pkg/seaice/seaice_preconditioner.F | 5 +- pkg/seaice/seaice_readparms.F | 2 +- .../1D_ocean_ice_column/code/SEAICE_OPTIONS.h | 209 ++++++++++++----- .../code_ad/SEAICE_OPTIONS.h | 212 +++++++++++++----- .../cpl_aim+ocn/code_ocn/SEAICE_OPTIONS.h | 191 +++++++++++----- .../code/SEAICE_OPTIONS.h | 205 ++++++++++++----- .../code_ad/SEAICE_OPTIONS.h | 205 ++++++++++++----- .../code_tap/SEAICE_OPTIONS.h | 207 ++++++++++++----- verification/lab_sea/code/SEAICE_OPTIONS.h | 205 ++++++++++++----- verification/lab_sea/code_ad/SEAICE_OPTIONS.h | 200 ++++++++++++----- .../lab_sea/code_tap/SEAICE_OPTIONS.h | 209 ++++++++++++----- .../offline_exf_seaice/code/SEAICE_OPTIONS.h | 210 ++++++++++++----- .../code_ad/SEAICE_OPTIONS.h | 205 ++++++++++++----- verification/seaice_itd/code/SEAICE_OPTIONS.h | 205 ++++++++++++----- .../seaice_obcs/code/SEAICE_OPTIONS.h | 205 ++++++++++++----- 50 files changed, 2367 insertions(+), 1145 deletions(-) diff --git a/doc/tag-index b/doc/tag-index index c93709d7e1..be6cc91462 100644 --- a/doc/tag-index +++ b/doc/tag-index @@ -1,6 +1,12 @@ Notes on tags used in MITgcmUV ============================== +o pkg/seaice: + - new option SEAICE_BGRID_DYNAMICS allows to remove/replace + SEAICE_ALLOW_DYNAMICS by ( SEAICE_CGRID or SEAICE_BGRID_DYNAMICS ) ; + - exclude more code if SEAICE_CGRID is undefined and clean up/rearrange + SEAICE.h to form more consistent groups and common-blocks ; + - update SEAICE_OPTIONS.h in verification experiments. o pkg/mdsio: - Add 2 new arguments (kLo & kHi) to pkg/mdsio Slice I/O routines that allow to select which levels in Slice-array to read-in or to write-out to file. diff --git a/model/src/the_main_loop.F b/model/src/the_main_loop.F index 151f4175aa..a4caf30de4 100644 --- a/model/src/the_main_loop.F +++ b/model/src/the_main_loop.F @@ -512,7 +512,7 @@ SUBROUTINE THE_MAIN_LOOP( myTime, myIter, myThid ) # endif /* ALLOW_BULKFORMULAE */ c-- # ifdef ALLOW_SEAICE -# ifdef SEAICE_ALLOW_DYNAMICS +# ifdef SEAICE_CGRID CADJ INIT comlev1_dynsol = COMMON,nchklev_1*MPSEUDOTIMESTEPS # ifdef SEAICE_LSR_ADJOINT_ITER CADJ INIT comlev1_lsr = COMMON,nchklev_1*MPSEUDOTIMESTEPS*SOLV_MAX_FIXED diff --git a/pkg/autodiff/autodiff_ini_model_io.F b/pkg/autodiff/autodiff_ini_model_io.F index 3940565913..161459818e 100644 --- a/pkg/autodiff/autodiff_ini_model_io.F +++ b/pkg/autodiff/autodiff_ini_model_io.F @@ -385,7 +385,7 @@ SUBROUTINE AUTODIFF_INI_MODEL_IO( myThid ) CALL MNC_CW_ADD_VATTR_TEXT('adhsnow', & 'coordinates','XC YC RC iter', myThid) c -# ifdef SEAICE_ALLOW_DYNAMICS +# ifdef SEAICE_CGRID CALL MNC_CW_ADD_VNAME('aduice', 'U_xy_Hn__-__t', 3,4, myThid) CALL MNC_CW_ADD_VATTR_TEXT('aduice', & 'units','[cost]/[m/s]', myThid) diff --git a/pkg/seaice/SEAICE.h b/pkg/seaice/SEAICE.h index 79ef9b0616..1fa0ae6fe2 100644 --- a/pkg/seaice/SEAICE.h +++ b/pkg/seaice/SEAICE.h @@ -36,6 +36,13 @@ C SIMaskU/V :: land-sea mask at U/V-points (copies of maskW/S(k=kSrf)) _RL HEFFM (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL SIMaskU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL SIMaskV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + +#if ( defined SEAICE_CGRID || defined SEAICE_BGRID_DYNAMICS ) + COMMON/ARRAYMETRIC/ k1AtC, k2AtC + _RS k1AtC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RS k2AtC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) +#endif + #ifdef SEAICE_CGRID COMMON/ARRAYC/ seaiceMaskU, seaiceMaskV C dynamic masks (depend on area) @@ -43,32 +50,46 @@ C dynamic masks (depend on area) _RL seaiceMaskV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C k1/2AtZ :: coefficients at C and Z points C k1/2AtC for metric terms in U/V ice equations. - COMMON/ARRAYCMETRIC/ k1AtC, k1AtZ, k2AtC, k2AtZ - _RS k1AtC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + COMMON/ARRAYCMETRIC/ k1AtZ, k2AtZ _RS k1AtZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RS k2AtC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RS k2AtZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) -#else +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS C UVM :: B-grid velocity-point mask COMMON/ARRAYB/ UVM _RS UVM (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C k1/2AtC/U/V :: coefficients at C, U, and V points C for metric terms in U/V ice equations. - COMMON/ARRAYBMETRIC/ - & k1AtC, k1AtU, k1AtV, k2AtC, k2AtU, k2AtV - _RS k1AtC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + COMMON/ARRAYBMETRIC/ k1AtU, k1AtV, k2AtU, k2AtV _RS k1AtU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RS k1AtV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RS k2AtC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RS k2AtU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RS k2AtV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) -#endif /* SEAICE_CGRID */ +#endif /* SEAICE_BGRID_DYNAMICS */ C-- Dynamical variables - COMMON/SEAICE_DYNVARS_1/AREA,HEFF,HSNOW,UICE,VICE + COMMON/SEAICE_DYNVARS_1/ + & AREA, HEFF, HSNOW, UICE, VICE + _RL AREA (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL HEFF (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL HSNOW (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL UICE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL VICE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) +C DWATN :: (linear) ice-ocean drag coefficient +C ( units of [rho|u|] = kg/m^2/s ) +C u/vIceNm1 :: sea ice drift velocities of previous timestep (m/s) + COMMON/SEAICE_DYNVARS_2/ + & DWATN, + & uIceNm1, vIceNm1 + _RL DWATN (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL uIceNm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL vIceNm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + #ifdef SEAICE_ITD - & ,AREAITD,HEFFITD,HSNOWITD, - & opnWtrFrac, fw2ObyRidge + COMMON/SEAICE_DYNVARS_ITD/ + & AREAITD, HEFFITD, HSNOWITD, + & opnWtrFrac, fw2ObyRidge C Fields for dynamic ice thickness distribution (ITD) C AREAITD :: area classes C HEFFITD :: ice thickenss classes (in meters) @@ -82,12 +103,28 @@ C ocean during ridging _RL opnWtrFrac (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL fw2ObyRidge(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) #endif /* SEAICE_ITD */ - _RL AREA (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL HEFF (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL HSNOW (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL UICE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL VICE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) +#ifdef SEAICE_CGRID +C stressDivergenceX/Y :: div of (vert. integr.) stress tensor (N/m^2) + COMMON /SEAICE_STRESSDIV/ + & stressDivergenceX, stressDivergenceY + _RL stressDivergenceX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL stressDivergenceY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) +# ifdef SEAICE_ALLOW_EVP +C-- Additional fields needed by the EVP solver: +C (vertically integrated) stress tensor, with diagonal terms sigma11/22 +C seaice_sigma1 :: sigma11+sigma22, defined at C-points (N/m) +C seaice_sigma2 :: sigma11-sigma22, defined at C-points (N/m) +C seaice_sigma12 :: off-diagonal term, defined at Z-points (N/m) + COMMON /SEAICE_EVP_FIELDS/ + & seaice_sigma1, seaice_sigma2, seaice_sigma12 + _RL seaice_sigma1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL seaice_sigma2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL seaice_sigma12 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) +# endif /* SEAICE_ALLOW_EVP */ +#endif + +#if ( defined SEAICE_CGRID || defined SEAICE_BGRID_DYNAMICS ) C ETA, etaZ :: shear viscosity as C-points, at Z-points (N s/m = kg/s) C ZETA, zetaA :: bulk viscosity at C-points, at Z-points C PRESS :: maximum vertically integrated ice strength/pressure (N/m) @@ -96,46 +133,39 @@ C deltaC :: deformation rate tensor invariant, for VP sea ice C = sqrt( (e11+e22)**2 + (1/e)*(e11-e22)**2 + 4*e12**2) ) C FORCEX/Y :: momentum forcing C ( units of [rho * h * u / deltaT] = kg/m/s^2 ) -C u/vIceNm1 :: sea ice drift velocities of previous timestep (m/s) +C tensileStrFac :: factor k to compute the maximal tensile stress k*PRESS0 COMMON/SEAICE_DYNVARS_3/ - & ETA,etaZ,ZETA,zetaZ,PRESS, e11, e22, e12, deltaC, - & FORCEX,FORCEY, - & uIceNm1, vIceNm1 + & ETA, etaZ, ZETA, zetaZ, PRESS, tensileStrFac, + & e11, e22, e12, deltaC, + & FORCEX,FORCEY + _RL ETA (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL etaZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL ZETA (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL zetaZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C ice strength/pressure term _RL PRESS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL tensileStrFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C strain rate tensor _RL e11 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL e22 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL e12 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) -C deformation rate tensor invariant, for viscous plastic sea ice = -C sqrt[(e11**2+e22**2)*(1+1/e**2)+ 4./e**2*e12C**2 + 2*e11*e22*(1-1/e**2)) _RL deltaC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C _RL FORCEX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL FORCEY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL uIceNm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL vIceNm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) -C DWATN :: (linear) ice-ocean drag coefficient -C ( units of [rho|u|] = kg/m^2/s ) C PRESS0 :: maximal compressive stress/strength (N/m) C FORCEX/Y0 :: external momentum forcing fields (part of FORCEX/Y) C SEAICE_zMax/zMin :: maximum/minimum bulk viscosities -C tensileStrFac :: factor k to compute the maximal tensile stress k*PRESS0 COMMON/SEAICE_DYNVARS_4/ - & DWATN, PRESS0, FORCEX0, FORCEY0, - & SEAICE_zMax, SEAICE_zMin, tensileStrFac - _RL DWATN (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL PRESS0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL FORCEX0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL FORCEY0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL SEAICE_zMax (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL SEAICE_zMin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL tensileStrFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + & PRESS0, FORCEX0, FORCEY0, SEAICE_zMax, SEAICE_zMin + _RL PRESS0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL FORCEX0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL FORCEY0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL SEAICE_zMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL SEAICE_zMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) +#endif #ifdef SEAICE_CGRID C seaiceMassC/U/V :: mass (ice+snow) at C/U/V-points ( kg/m^2 ) @@ -144,11 +174,6 @@ C seaiceMassC/U/V :: mass (ice+snow) at C/U/V-points ( kg/m^2 ) _RL seaiceMassC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL seaiceMassU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL seaiceMassV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) -C stressDivergenceX/Y :: div of (vert. integr.) stress tensor (N/m^2) - COMMON /SEAICE_STRESSDIV/ - & stressDivergenceX, stressDivergenceY - _RL stressDivergenceX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL stressDivergenceY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) # ifdef SEAICE_ALLOW_FREEDRIFT C u/vice_fd :: free drift velocities (m/s) COMMON /SEAICE_FD_FIELDS/ @@ -157,26 +182,13 @@ C u/vice_fd :: free drift velocities (m/s) _RL vice_fd (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) # endif -# ifdef SEAICE_ALLOW_EVP -C-- Additional fields needed by the EVP solver: -C (vertically integrated) stress tensor, with diagonal terms sigma11/22 -C seaice_sigma1 :: sigma11+sigma22, defined at C-points (N/m) -C seaice_sigma2 :: sigma11-sigma22, defined at C-points (N/m) -C seaice_sigma12 :: off-diagonal term, defined at Z-points (N/m) - COMMON /SEAICE_EVP_FIELDS/ - & seaice_sigma1, seaice_sigma2, seaice_sigma12 - _RL seaice_sigma1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL seaice_sigma2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL seaice_sigma12 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) -# endif /* SEAICE_ALLOW_EVP */ - # ifdef SEAICE_ALLOW_BOTTOMDRAG C CbobC :: (linear) bottom drag coefficient for basals stress param. COMMON/SEAICE_BOTTOMDRAG/ CbotC _RL CbotC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) # endif /* SEAICE_ALLOW_BOTTOMDRAG */ -# if (defined SEAICE_ALLOW_JFNK) || (defined SEAICE_ALLOW_KRYLOV) +# if ( defined SEAICE_ALLOW_JFNK || defined SEAICE_ALLOW_KRYLOV ) C diagnostics for the JFNK and Krylov solver INTEGER totalNewtonIters INTEGER totalNewtonFails @@ -187,22 +199,26 @@ C diagnostics for the JFNK and Krylov solver & totalNewtonIters, totalNewtonFails, & totalKrylovIters, totalKrylovFails, & totalJFNKtimeSteps +C Scalar product used in FGMRES needs a metric INTEGER nVec PARAMETER ( nVec=2*sNx*sNy ) _RL scalarProductMetric( nVec, 1, nSx, nSy ) COMMON /SEAICE_KRYLOV_RL/ scalarProductMetric # endif /* SEAICE_ALLOW_JFNK or SEAICE_ALLOW_KRYLOV */ -#else /* ndef SEAICE_CGRID */ +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS C AMASS :: sea ice mass C DAIRN :: (linear) atmosphere-ice drag coefficient -C uIceC :: average of UICE between last two time steps -C vIceC :: average of VICE between last two time steps - COMMON/SEAICE_DYNVARS_BGRID/ AMASS, DAIRN, uIceC, vIceC +C u/vIceC has been renamed to u/vIceC to avoid conficts with C-grid code +C uIceB :: average of UICE between last two time steps +C vIceB :: average of VICE between last two time steps + COMMON/SEAICE_DYNVARS_BGRID/ AMASS, DAIRN, uIceB, vIceB _RL AMASS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL DAIRN (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL uIceC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) - _RL vIceC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL uIceB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL vIceB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) COMMON/WIND_STRESS_OCE/WINDX,WINDY C WINDX - zonal wind stress over water at C points @@ -218,7 +234,7 @@ C GWATX/Y :: geostrophic ocean velocities C-- KGEO Level used as a proxy for geostrophic velocity. COMMON/SEAICE_KGEO/KGEO INTEGER KGEO (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) -#endif /* SEAICE_CGRID */ +#endif /* SEAICE_BGRID_DYNAMICS */ COMMON/SEAICE_REG_NEG/d_HEFFbyNEG,d_HSNWbyNEG C The change of mean ice thickness due to out-of-bounds values following diff --git a/pkg/seaice/SEAICE_OPTIONS.h b/pkg/seaice/SEAICE_OPTIONS.h index 425e9afcff..b9f92251f8 100644 --- a/pkg/seaice/SEAICE_OPTIONS.h +++ b/pkg/seaice/SEAICE_OPTIONS.h @@ -79,10 +79,19 @@ C computed the "call diagnostics_fill" statement is commented out. # undef ALLOW_SITRACER_DEBUG_DIAG #endif /* ALLOW_SITRACER */ +C-- Allow sea-ice dynamic code. These options are 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). + 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. +C discretization should still work but it is not longer actively +C tested and supported. Define this flag to compile it. It cannot be +C defined together with SEAICE_CGRID +#undef SEAICE_BGRID_DYNAMICS + +C-- The following flag should always be set in order to use C the +C-- operational C-grid discretization. #define SEAICE_CGRID #ifdef SEAICE_CGRID @@ -174,7 +183,9 @@ 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 */ +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS C-- Options for the B-grid version only: C- By default for B-grid dynamics solver wind stress under sea-ice is @@ -192,7 +203,7 @@ 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 */ +#endif /* SEAICE_BGRID_DYNAMICS */ C-- Some regularisations C- When set limit the Ice-Loading to mass of 1/5 of Surface ocean grid-box @@ -242,12 +253,6 @@ 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 diff --git a/pkg/seaice/dynsolver.F b/pkg/seaice/dynsolver.F index 0a4b92b5f5..49e12788f6 100644 --- a/pkg/seaice/dynsolver.F +++ b/pkg/seaice/dynsolver.F @@ -10,8 +10,8 @@ SUBROUTINE DYNSOLVER( myTime, myIter, myThid ) C *==========================================================* C | SUBROUTINE DYNSOLVER | -C | o Ice dynamics using LSR solver | -C | Zhang and Hibler, JGR, 102, 8691-8702, 1997 | +C | o B-grid version of ice dynamics using LSR solver | +C | Zhang and Hibler, JGR, 102, 8691-8702, 1997 | C *==========================================================* C *==========================================================* IMPLICIT NONE @@ -42,7 +42,7 @@ SUBROUTINE DYNSOLVER( myTime, myIter, myThid ) INTEGER myThid CEndOfInterface -#ifndef SEAICE_CGRID +#ifdef SEAICE_BGRID_DYNAMICS #ifdef EXPLICIT_SSH_SLOPE #include "SURFACE.h" @@ -233,8 +233,6 @@ SUBROUTINE DYNSOLVER( myTime, myIter, myThid ) ENDDO ENDDO -#ifdef SEAICE_ALLOW_DYNAMICS - IF ( SEAICEuseDYNAMICS ) THEN #ifdef ALLOW_AUTODIFF_TAMC @@ -255,8 +253,8 @@ SUBROUTINE DYNSOLVER( myTime, myIter, myThid ) DO i=1-OLx,sNx+OLx UICENM1(i,j,bi,bj)=UICE(i,j,bi,bj) VICENM1(i,j,bi,bj)=VICE(i,j,bi,bj) - UICEC(i,j,bi,bj)=UICE(i,j,bi,bj) - VICEC(i,j,bi,bj)=VICE(i,j,bi,bj) + UICEB(i,j,bi,bj)=UICE(i,j,bi,bj) + VICEB(i,j,bi,bj)=VICE(i,j,bi,bj) ENDDO ENDDO ENDDO @@ -276,8 +274,8 @@ SUBROUTINE DYNSOLVER( myTime, myIter, myThid ) DO i=1-OLx,sNx+OLx UICE(i,j,bi,bj)=HALF*(UICE(i,j,bi,bj)+UICENM1(i,j,bi,bj)) VICE(i,j,bi,bj)=HALF*(VICE(i,j,bi,bj)+VICENM1(i,j,bi,bj)) - UICEC(i,j,bi,bj)=UICE(i,j,bi,bj) - VICEC(i,j,bi,bj)=VICE(i,j,bi,bj) + UICEB(i,j,bi,bj)=UICE(i,j,bi,bj) + VICEB(i,j,bi,bj)=VICE(i,j,bi,bj) ENDDO ENDDO ENDDO @@ -289,7 +287,6 @@ SUBROUTINE DYNSOLVER( myTime, myIter, myThid ) cdm c$taf store uice,vice = comlev1, key=ikey_dynamics ENDIF -#endif /* SEAICE_ALLOW_DYNAMICS */ #ifdef ALLOW_DIAGNOSTICS IF ( useDiagnostics ) THEN @@ -302,7 +299,6 @@ SUBROUTINE DYNSOLVER( myTime, myIter, myThid ) C Calculate ocean surface stress CALL OSTRES ( COR_ICE, myThid ) -#ifdef SEAICE_ALLOW_DYNAMICS #ifdef SEAICE_ALLOW_CLIPVELS IF ( SEAICEuseDYNAMICS .AND. SEAICE_clipVelocities) THEN #ifdef ALLOW_AUTODIFF_TAMC @@ -347,8 +343,7 @@ SUBROUTINE DYNSOLVER( myTime, myIter, myThid ) ENDIF #endif /* SEAICE_ALLOW_CLIPVELS */ -#endif /* SEAICE_ALLOW_DYNAMICS */ -#endif /* NOT SEAICE_CGRID */ +#endif /* SEAICE_BGRID_DYNAMICS */ RETURN END diff --git a/pkg/seaice/lsr.F b/pkg/seaice/lsr.F index 886788f2ce..dbbee44679 100644 --- a/pkg/seaice/lsr.F +++ b/pkg/seaice/lsr.F @@ -40,8 +40,7 @@ SUBROUTINE lsr( ilcall, myThid ) INTEGER myThid CEndOfInterface -#ifndef SEAICE_CGRID -#ifdef SEAICE_ALLOW_DYNAMICS +#ifdef SEAICE_BGRID_DYNAMICS C === Local variables === C i,j,bi,bj - Loop counters @@ -168,30 +167,30 @@ SUBROUTINE lsr( ilcall, myThid ) DO j=1-OLy,sNy+OLy-1 DO i=1-OLx,sNx+OLx-1 ucLoc(I,J,bi,bj) = 0.25 _d 0 * ( - & + uIceC(I+1,J,bi,bj) + uIceC(I+1,J+1,bi,bj) - & + uIceC(I ,J,bi,bj) + uIceC(I, J+1,bi,bj) ) + & + uIceB(I+1,J,bi,bj) + uIceB(I+1,J+1,bi,bj) + & + uIceB(I ,J,bi,bj) + uIceB(I, J+1,bi,bj) ) vcLoc(I,J,bi,bj) = 0.25 _d 0 * ( - & + vIceC(I+1,J,bi,bj) + vIceC(I+1,J+1,bi,bj) - & + vIceC(I ,J,bi,bj) + vIceC(I, J+1,bi,bj) ) + & + vIceB(I+1,J,bi,bj) + vIceB(I+1,J+1,bi,bj) + & + vIceB(I ,J,bi,bj) + vIceB(I, J+1,bi,bj) ) ENDDO ENDDO DO J=1-OLy,sNy+OLy-1 DO I=1-OLy,sNx+OLy-1 e11loc = 0.5 _d 0 * _recip_dxF(I,J,bi,bj) * - & (uIceC(I+1,J+1,bi,bj)+uIceC(I+1,J,bi,bj) - & -uIceC(I, J+1,bi,bj)-uIceC(I, J,bi,bj)) + & (uIceB(I+1,J+1,bi,bj)+uIceB(I+1,J,bi,bj) + & -uIceB(I, J+1,bi,bj)-uIceB(I, J,bi,bj)) & + vcLoc(I,J,bi,bj) * k2AtC(I,J,bi,bj) e22loc = 0.5 _d 0 * _recip_dyF(I,J,bi,bj) * - & (vIceC(I+1,J+1,bi,bj)+vIceC(I,J+1,bi,bj) - & -vIceC(I+1,J, bi,bj)-vIceC(I,J, bi,bj)) + & (vIceB(I+1,J+1,bi,bj)+vIceB(I,J+1,bi,bj) + & -vIceB(I+1,J, bi,bj)-vIceB(I,J, bi,bj)) & + ucLoc(I,J,bi,bj) * k1AtC(I,J,bi,bj) e12loc = 0.5 _d 0*( & 0.5 _d 0 * _recip_dyF(I,J,bi,bj) * - & (uIceC(I+1,J+1,bi,bj)+uIceC(I,J+1,bi,bj) - & -uIceC(I+1,J, bi,bj)-uIceC(I,J, bi,bj)) + & (uIceB(I+1,J+1,bi,bj)+uIceB(I,J+1,bi,bj) + & -uIceB(I+1,J, bi,bj)-uIceB(I,J, bi,bj)) & + 0.5 _d 0 * _recip_dxF(I,J,bi,bj) * - & (vIceC(I+1,J+1,bi,bj)+vIceC(I+1,J,bi,bj) - & -vIceC(I, J+1,bi,bj)-vIceC(I, J,bi,bj)) + & (vIceB(I+1,J+1,bi,bj)+vIceB(I+1,J,bi,bj) + & -vIceB(I, J+1,bi,bj)-vIceB(I, J,bi,bj)) & - vcLoc(I,J,bi,bj)*k1AtC(I,J,bi,bj) & - ucLoc(I,J,bi,bj)*k2AtC(I,J,bi,bj) ) C NOW EVALUATE VISCOSITIES @@ -412,13 +411,13 @@ SUBROUTINE lsr( ilcall, myThid ) & * (vcLoc(I,J,bi,bj)-vcLoc(I,J-1,bi,bj)) & * _recip_dyC(I,J,bi,bj) & + (zetaV(I,J)+etaV(I,J))*k2atV(I,J,bi,bj) - & * 0.5 _d 0 * (vIceC(I,J,bi,bj)+vIceC(I+1,J,bi,bj)) + & * 0.5 _d 0 * (vIceB(I,J,bi,bj)+vIceB(I+1,J,bi,bj)) & - 0.5 _d 0 * (PRESS(I,J,bi,bj)+PRESS(I,J-1,bi,bj)) sig12(I,J) = etaV(I,J) & * (ucLoc(I,J,bi,bj)-ucLoc(I,J-1,bi,bj)) & * _recip_dyC(I,J,bi,bj) & - etaV(I,J) * k2AtV(I,J,bi,bj) - & * 0.5 _d 0 * (uIceC(I,J,bi,bj)+uIceC(I+1,J,bi,bj)) + & * 0.5 _d 0 * (uIceB(I,J,bi,bj)+uIceB(I+1,J,bi,bj)) ENDDO ENDDO DO j=0,sNy @@ -429,25 +428,25 @@ SUBROUTINE lsr( ilcall, myThid ) & * (ucLoc(I,J,bi,bj)-ucLoc(I-1,J,bi,bj)) & * _recip_dxC(I,J,bi,bj) & + (zetaU(I,J)+etaU(I,J))*k2atU(I,J,bi,bj) - & * 0.5 _d 0 * (uIceC(I,J,bi,bj)+uIceC(I,J+1,bi,bj)) + & * 0.5 _d 0 * (uIceB(I,J,bi,bj)+uIceB(I,J+1,bi,bj)) & - 0.5 _d 0 * (PRESS(I,J,bi,bj)+PRESS(I-1,J,bi,bj)) sig21(I,J) = etaU(I,J) & * (vcLoc(I,J,bi,bj)-vcLoc(I-1,J,bi,bj)) & * _recip_dxC(I,J,bi,bj) & - etaU(I,J) * k1AtU(I,J,bi,bj) - & * 0.5 _d 0 * (vIceC(I,J,bi,bj)+vIceC(I,J+1,bi,bj)) + & * 0.5 _d 0 * (vIceB(I,J,bi,bj)+vIceB(I,J+1,bi,bj)) ENDDO ENDDO DO j=1,sNy DO i=1,sNx - rhsU(I,J,bi,bj) = DRAGA(I,J,bi,bj)*vIceC(I,J,bi,bj) + rhsU(I,J,bi,bj) = DRAGA(I,J,bi,bj)*vIceB(I,J,bi,bj) & +FORCEX(I,J,bi,bj) & + ( _dyC(I, J, bi,bj) * sig11(I, J ) & - _dyC(I-1,J, bi,bj) * sig11(I-1,J ) & + _dxC(I, J, bi,bj) * sig21(I, J ) & - _dxC(I, J-1,bi,bj) * sig21(I, J-1) ) & * recip_rAz(I,J,bi,bj) * UVM(I,J,bi,bj) - rhsV(I,J,bi,bj) = - DRAGA(I,J,bi,bj)*uIceC(I,J,bi,bj) + rhsV(I,J,bi,bj) = - DRAGA(I,J,bi,bj)*uIceB(I,J,bi,bj) & +FORCEY(I,J,bi,bj) & + ( _dyC(I, J, bi,bj) * sig12(I, J ) & - _dyC(I-1,J, bi,bj) * sig12(I-1,J ) @@ -746,14 +745,14 @@ SUBROUTINE lsr( ilcall, myThid ) DO J=1-OLy,sNy+OLy-1 DO I=1-OLx,sNx+OLx-1 dVdy(I,J) = 0.5 _d 0 * ( - & ( VICEC(I+1,J+1,bi,bj) - VICEC(I+1,J ,bi,bj) ) + & ( vIceB(I+1,J+1,bi,bj) - vIceB(I+1,J ,bi,bj) ) & * _recip_dyG(I+1,J,bi,bj) - & +(VICEC(I ,J+1,bi,bj) - VICEC(I ,J ,bi,bj) ) + & +(vIceB(I ,J+1,bi,bj) - vIceB(I ,J ,bi,bj) ) & * _recip_dyG(I, J,bi,bj) ) dVdx(I,J) = 0.5 _d 0 * ( - & ( VICEC(I+1,J+1,bi,bj) - VICEC(I ,J+1,bi,bj) ) + & ( vIceB(I+1,J+1,bi,bj) - vIceB(I ,J+1,bi,bj) ) & * _recip_dxG(I,J+1,bi,bj) - & +(VICEC(I+1,J ,bi,bj) - VICEC(I ,J ,bi,bj) ) + & +(vIceB(I+1,J ,bi,bj) - vIceB(I ,J ,bi,bj) ) & * _recip_dxG(I,J, bi,bj) ) ENDDO ENDDO @@ -761,15 +760,15 @@ SUBROUTINE lsr( ilcall, myThid ) DO j=1-OLy,sNy+OLy-1 DO i=1-OLx,sNx+OLx-1 vz(i,j) = quart * ( - & vicec(i,j,bi,bj) + vicec(i+1,j,bi,bj) ) + & vIceB(i,j,bi,bj) + vIceB(i+1,j,bi,bj) ) vz(i,j)= vz(i,j) + quart * ( - & vicec(i,j+1,bi,bj) + vicec(i+1,j+1,bi,bj) ) + & vIceB(i,j+1,bi,bj) + vIceB(i+1,j+1,bi,bj) ) ENDDO ENDDO #endif DO J=1,sNy DO I=1,sNx - rhsU(I,J,bi,bj)=DRAGA(I,J,bi,bj)*VICEC(I,J,bi,bj) + rhsU(I,J,bi,bj)=DRAGA(I,J,bi,bj)*vIceB(I,J,bi,bj) & +FORCEX(I,J,bi,bj) #ifdef SEAICE_TEST & + ( 0.5 _d 0 * @@ -795,18 +794,18 @@ SUBROUTINE lsr( ilcall, myThid ) & & -(etaPlusZeta(I ,J ,bi,bj)+etaPlusZeta(I ,J-1,bi,bj) & -etaPlusZeta(I-1,J-1,bi,bj)-etaPlusZeta(I-1,J ,bi,bj)) - & * VICEC(I,J,bi,bj) + & * vIceB(I,J,bi,bj) & * _tanPhiAtV(I,J,bi,bj) & * 0.5 _d 0 * _recip_dxV(I,J,bi,bj)*recip_rSphere & & -(ETAMEAN(I,J,bi,bj)+ZETAMEAN(I,J,bi,bj)) - & *(VICEC(I+1,J,bi,bj) - VICEC(I-1,J,bi,bj)) + & *(vIceB(I+1,J,bi,bj) - vIceB(I-1,J,bi,bj)) & * _tanPhiAtV(I,J,bi,bj) & * 1.0 _d 0 /( _dxG(I,J,bi,bj) + _dxG(I-1,J,bi,bj) ) & *recip_rSphere & & -ETAMEAN(I,J,bi,bj) - & *(VICEC(I+1,J,bi,bj) - VICEC(I-1,J,bi,bj)) + & *(vIceB(I+1,J,bi,bj) - vIceB(I-1,J,bi,bj)) & *TWO* _tanPhiAtV(I,J,bi,bj) & * 1.0 _d 0 /( _dxG(I,J,bi,bj) + _dxG(I-1,J,bi,bj) ) & *recip_rSphere @@ -989,14 +988,14 @@ SUBROUTINE lsr( ilcall, myThid ) DO J=1-OLy,sNy+OLy-1 DO I=1-OLx,sNx+OLx-1 dUdx(I,J) = 0.5 _d 0 * ( - & ( UICEC(I+1,J+1,bi,bj) - UICEC(I ,J+1,bi,bj) ) + & ( uIceB(I+1,J+1,bi,bj) - uIceB(I ,J+1,bi,bj) ) & * _recip_dxG(I,J+1,bi,bj) - & +(UICEC(I+1,J ,bi,bj) - UICEC(I ,J ,bi,bj) ) + & +(uIceB(I+1,J ,bi,bj) - uIceB(I ,J ,bi,bj) ) & * _recip_dxG(I,J ,bi,bj) ) dUdy(I,J) = 0.5 _d 0 * ( - & ( UICEC(I+1,J+1,bi,bj) - UICEC(I+1,J ,bi,bj) ) + & ( uIceB(I+1,J+1,bi,bj) - uIceB(I+1,J ,bi,bj) ) & * _recip_dyG(I+1,J,bi,bj) - & +(UICEC(I ,J+1,bi,bj) - UICEC(I ,J ,bi,bj) ) + & +(uIceB(I ,J+1,bi,bj) - uIceB(I ,J ,bi,bj) ) & * _recip_dyG(I, J,bi,bj) ) ENDDO ENDDO @@ -1004,15 +1003,15 @@ SUBROUTINE lsr( ilcall, myThid ) DO j=1-OLy,sNy+OLy-1 DO i=1-OLx,sNx+OLx-1 uz(i,j) = quart * ( - & uicec(i,j,bi,bj) + uicec(i+1,j,bi,bj) ) + & uIceB(i,j,bi,bj) + uIceB(i+1,j,bi,bj) ) uz(i,j)= uz(i,j) + quart * ( - & uicec(i,j+1,bi,bj) + uicec(i+1,j+1,bi,bj) ) + & uIceB(i,j+1,bi,bj) + uIceB(i+1,j+1,bi,bj) ) ENDDO ENDDO #endif DO J=1,sNy DO I=1,sNx - rhsV(I,J,bi,bj)=-DRAGA(I,J,bi,bj)*UICEC(I,J,bi,bj) + rhsV(I,J,bi,bj)=-DRAGA(I,J,bi,bj)*uIceB(I,J,bi,bj) & +FORCEY(I,J,bi,bj) & #ifdef SEAICE_TEST @@ -1039,15 +1038,15 @@ SUBROUTINE lsr( ilcall, myThid ) & & +(ETA(I ,J ,bi,bj) + ETA(I ,J-1,bi,bj) & -ETA(I-1,J-1,bi,bj) - ETA(I-1,J ,bi,bj)) - & * UICEC(I,J,bi,bj) + & * uIceB(I,J,bi,bj) & * _tanPhiAtV(I,J,bi,bj) & * 0.5 _d 0 * _recip_dxV(I,J,bi,bj)*recip_rSphere & +ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj) - & *(UICEC(I+1,J,bi,bj)-UICEC(I-1,J,bi,bj)) + & *(uIceB(I+1,J,bi,bj)-uIceB(I-1,J,bi,bj)) & * 0.5 _d 0 * _recip_dxV(I,J,bi,bj)*recip_rSphere & & +ETAMEAN(I,J,bi,bj)*TWO * _tanPhiAtV(I,J,bi,bj) - & *(UICEC(I+1,J,bi,bj)-UICEC(I-1,J,bi,bj)) + & *(uIceB(I+1,J,bi,bj)-uIceB(I-1,J,bi,bj)) & * 1. _d 0 /( _dxG(I,J,bi,bj) + _dxG(I-1,J,bi,bj)) & *recip_rSphere VRT1(I,J,bi,bj)= 0.5 _d 0 * ( @@ -1197,8 +1196,7 @@ SUBROUTINE lsr( ilcall, myThid ) ENDDO ENDDO -#endif /* SEAICE_ALLOW_DYNAMICS */ -#endif /* SEAICE_CGRID */ +#endif /* SEAICE_BGRID_DYNAMICS */ RETURN END diff --git a/pkg/seaice/ostres.F b/pkg/seaice/ostres.F index e9f012c3fb..3692635e55 100644 --- a/pkg/seaice/ostres.F +++ b/pkg/seaice/ostres.F @@ -23,15 +23,15 @@ SUBROUTINE ostres( COR_ICE, myThid ) INTEGER myThid CEndOfInterface -#ifndef SEAICE_CGRID +#ifdef SEAICE_BGRID_DYNAMICS C === Local variables === C i,j,bi,bj - Loop counters INTEGER i, j, bi, bj _RL SINWIN, COSWIN, SINWAT, COSWAT -#ifdef SEAICE_BICE_STRESS +# ifdef SEAICE_BICE_STRESS _RL fuIce, fvIce -#endif +# endif c introduce turning angle (default is zero) SINWIN=SIN(SEAICE_airTurnAngle*deg2rad) @@ -42,7 +42,7 @@ SUBROUTINE ostres( COR_ICE, myThid ) C-- Update overlap regions CALL EXCH_UV_XY_RL(WINDX, WINDY, .TRUE., myThid) -#ifndef SEAICE_EXTERNAL_FLUXES +# ifndef SEAICE_EXTERNAL_FLUXES C-- Interpolate wind stress (N/m^2) from South-West B-grid C to South-West C-grid for forcing ocean model. DO bj=myByLo(myThid),myByHi(myThid) @@ -58,9 +58,9 @@ SUBROUTINE ostres( COR_ICE, myThid ) ENDDO ENDDO CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid) -#endif /* ifndef SEAICE_EXTERNAL_FLUXES */ +# endif /* ifndef SEAICE_EXTERNAL_FLUXES */ -#ifdef SEAICE_BICE_STRESS +# ifdef SEAICE_BICE_STRESS C-- Compute ice-affected wind stress DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) @@ -91,8 +91,8 @@ SUBROUTINE ostres( COR_ICE, myThid ) ENDDO ENDDO CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid) -#endif /* SEAICE_BICE_STRESS */ -#endif /* not SEAICE_CGRID */ +# endif /* SEAICE_BICE_STRESS */ +#endif /* SEAICE_BCGRID_DYNAMICS */ RETURN END diff --git a/pkg/seaice/seaice_ad_check_lev1_dir.h b/pkg/seaice/seaice_ad_check_lev1_dir.h index 9dc3a8f78f..ce280d1192 100644 --- a/pkg/seaice/seaice_ad_check_lev1_dir.h +++ b/pkg/seaice/seaice_ad_check_lev1_dir.h @@ -16,7 +16,7 @@ CADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE totphihyd = comlev1, key=ikey_dynamics, kind=isbyte cphCADJ STORE runoff = comlev1, key=ikey_dynamics, kind=isbyte -# ifdef SEAICE_ALLOW_DYNAMICS +# ifdef SEAICE_CGRID cphCADJ STORE zeta = comlev1, key=ikey_dynamics, kind=isbyte cphCADJ STORE eta = comlev1, key=ikey_dynamics, kind=isbyte # endif @@ -25,17 +25,13 @@ cphCADJ STORE eta = comlev1, key=ikey_dynamics, kind=isbyte #ifdef SEAICE_CGRID CADJ STORE stressdivergencex = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE stressdivergencey = comlev1, key=ikey_dynamics, kind=isbyte -#endif -#ifdef SEAICE_ALLOW_DYNAMICS -# ifdef SEAICE_CGRID CADJ STORE etan = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE dwatn = comlev1, key=ikey_dynamics, kind=isbyte -# ifdef SEAICE_ALLOW_BOTTOMDRAG +# ifdef SEAICE_ALLOW_BOTTOMDRAG CADJ STORE cbotc = comlev1, key=ikey_dynamics, kind=isbyte -# endif +# endif CADJ STORE seaicemasku = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE seaicemaskv = comlev1, key=ikey_dynamics, kind=isbyte -# endif # ifdef SEAICE_ALLOW_EVP CADJ STORE seaice_sigma1 = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE seaice_sigma2 = comlev1, key=ikey_dynamics, kind=isbyte diff --git a/pkg/seaice/seaice_ad_check_lev2_dir.h b/pkg/seaice/seaice_ad_check_lev2_dir.h index 8b7c03b319..f4ce8b9a8c 100644 --- a/pkg/seaice/seaice_ad_check_lev2_dir.h +++ b/pkg/seaice/seaice_ad_check_lev2_dir.h @@ -13,10 +13,10 @@ CADJ STORE HSNOW = tapelev2, key = ilev_2 CADJ STORE RUNOFF = tapelev2, key = ilev_2 CADJ STORE UICE = tapelev2, key = ilev_2 CADJ STORE VICE = tapelev2, key = ilev_2 -CADJ STORE ZETA = tapelev2, key = ilev_2 -CADJ STORE ETA = tapelev2, key = ilev_2 CADJ STORE TICES = tapelev2, key = ilev_2 # ifdef SEAICE_CGRID +CADJ STORE ZETA = tapelev2, key = ilev_2 +CADJ STORE ETA = tapelev2, key = ilev_2 CADJ STORE dwatn = tapelev2, key = ilev_2 # ifdef SEAICE_ALLOW_BOTTOMDRAG CADJ STORE cbotc = tapelev2, key = ilev_2 diff --git a/pkg/seaice/seaice_ad_check_lev3_dir.h b/pkg/seaice/seaice_ad_check_lev3_dir.h index ea4ff789de..31f5770192 100644 --- a/pkg/seaice/seaice_ad_check_lev3_dir.h +++ b/pkg/seaice/seaice_ad_check_lev3_dir.h @@ -13,10 +13,10 @@ CADJ STORE HSNOW = tapelev3, key = ilev_3 CADJ STORE RUNOFF = tapelev3, key = ilev_3 CADJ STORE UICE = tapelev3, key = ilev_3 CADJ STORE VICE = tapelev3, key = ilev_3 -CADJ STORE ZETA = tapelev3, key = ilev_3 -CADJ STORE ETA = tapelev3, key = ilev_3 CADJ STORE TICES = tapelev3, key = ilev_3 # ifdef SEAICE_CGRID +CADJ STORE ZETA = tapelev3, key = ilev_3 +CADJ STORE ETA = tapelev3, key = ilev_3 CADJ STORE dwatn = tapelev3, key = ilev_3 # ifdef SEAICE_ALLOW_BOTTOMDRAG CADJ STORE cbotc = tapelev3, key = ilev_3 diff --git a/pkg/seaice/seaice_ad_check_lev4_dir.h b/pkg/seaice/seaice_ad_check_lev4_dir.h index 55cd8a8fdb..c2989e2915 100644 --- a/pkg/seaice/seaice_ad_check_lev4_dir.h +++ b/pkg/seaice/seaice_ad_check_lev4_dir.h @@ -13,10 +13,10 @@ CADJ STORE HSNOW = tapelev4, key = ilev_4 CADJ STORE RUNOFF = tapelev4, key = ilev_4 CADJ STORE UICE = tapelev4, key = ilev_4 CADJ STORE VICE = tapelev4, key = ilev_4 -CADJ STORE ZETA = tapelev4, key = ilev_4 -CADJ STORE ETA = tapelev4, key = ilev_4 CADJ STORE TICES = tapelev4, key = ilev_4 # ifdef SEAICE_CGRID +CADJ STORE ZETA = tapelev4, key = ilev_4 +CADJ STORE ETA = tapelev4, key = ilev_4 CADJ STORE dwatn = tapelev4, key = ilev_4 # ifdef SEAICE_ALLOW_BOTTOMDRAG CADJ STORE cbotc = tapelev4, key = ilev_4 diff --git a/pkg/seaice/seaice_ad_dump.F b/pkg/seaice/seaice_ad_dump.F index 92be4bceae..610c17a354 100644 --- a/pkg/seaice/seaice_ad_dump.F +++ b/pkg/seaice/seaice_ad_dump.F @@ -140,14 +140,14 @@ SUBROUTINE SEAICE_AD_DUMP( myTime, myIter, myThid ) & 12, doDump, dumpAdRecSi, myTime, myIter, myThid ) CALL DUMP_ADJ_XY( dumRS, adhsnow, 'ADJhsnow', 'ADJhsnow.', & 12, doDump, dumpAdRecSi, myTime, myIter, myThid ) -# ifdef SEAICE_ALLOW_DYNAMICS +# if ( defined SEAICE_CGRID || defined SEAICE_BGRID_DYNAMICS ) C IF (SEAICEuseDYNAMICS) THEN CALL DUMP_ADJ_XY_UV( & dumRS, aduice, 'ADJuice ', 'ADJuice.', & dumRS, advice, 'ADJvice ', 'ADJvice.', & 34, doDump, dumpAdRecSi, myTime, myIter, myThid) C ENDIF -# endif /* SEAICE_ALLOW_DYNAMICS */ +# endif IF ( doDump ) THEN #ifdef ALLOW_MNC @@ -170,7 +170,7 @@ SUBROUTINE SEAICE_AD_DUMP( myTime, myIter, myThid ) CALL COPY_ADVAR_OUTP( dumRS, adhsnow, var2Du, 1 , 12, myThid ) CALL MNC_CW_RL_W('D','adseaice',0,0,'adhsnow', & var2Du, myThid) -# ifdef SEAICE_ALLOW_DYNAMICS +# if ( defined SEAICE_CGRID || defined SEAICE_BGRID_DYNAMICS ) C IF (SEAICEuseDYNAMICS) THEN CALL COPY_AD_UV_OUTP( dumRS, dumRS, aduice, advice, & var2Du, var2Dv, 1, 34, myThid ) @@ -187,7 +187,7 @@ SUBROUTINE SEAICE_AD_DUMP( myTime, myIter, myThid ) & adheff, myThid) CALL MNC_CW_RL_W('D','adseaice',0,0,'adhsnow', & adhsnow, myThid) -# ifdef SEAICE_ALLOW_DYNAMICS +# if ( defined SEAICE_CGRID || defined SEAICE_BGRID_DYNAMICS ) C IF (SEAICEuseDYNAMICS) THEN CALL MNC_CW_RL_W('D','adseaice',0,0,'aduice', & aduice, myThid) diff --git a/pkg/seaice/seaice_advdiff.F b/pkg/seaice/seaice_advdiff.F index e9a6987551..55b5b6188f 100644 --- a/pkg/seaice/seaice_advdiff.F +++ b/pkg/seaice/seaice_advdiff.F @@ -90,39 +90,41 @@ SUBROUTINE SEAICE_ADVDIFF( _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) +#if ( defined SEAICE_ITD || !defined ALLOW_GENERIC_ADVDIFF ) CHARACTER*(MAX_LEN_MBUF) msgBuf +#endif CEOP 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 +C-- compute cell areas used by all tracers DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) -#ifndef SEAICE_CGRID /* not SEAICE_CGRID = BGRID */ + 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 + +#ifdef SEAICE_BGRID_DYNAMICS C-- hack to ensure backward compatibility: C average B-grid seaice velocity to C-grid + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) 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 */ +#endif #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE uc = comlev1, key = ikey_dynamics, kind=isbyte diff --git a/pkg/seaice/seaice_calc_ice_strength.F b/pkg/seaice/seaice_calc_ice_strength.F index e05918ae6b..2d8bba7ec2 100644 --- a/pkg/seaice/seaice_calc_ice_strength.F +++ b/pkg/seaice/seaice_calc_ice_strength.F @@ -45,6 +45,7 @@ SUBROUTINE SEAICE_CALC_ICE_STRENGTH( INTEGER myThid CEOP +#if ( defined SEAICE_CGRID || defined SEAICE_BGRID_DYNAMICS ) C !LOCAL VARIABLES: ==================================================== C === Local variables === C i,j,k :: inner loop counters @@ -192,6 +193,7 @@ SUBROUTINE SEAICE_CALC_ICE_STRENGTH( CML ENDDO CML ENDIF CML#endif /* SEAICE_CGRID */ +#endif /* SEAICE_CGRID or SEAICE_BGRID_DYNAMICS */ RETURN END diff --git a/pkg/seaice/seaice_calc_lhs.F b/pkg/seaice/seaice_calc_lhs.F index ee9a33dc65..9eea1c16a9 100644 --- a/pkg/seaice/seaice_calc_lhs.F +++ b/pkg/seaice/seaice_calc_lhs.F @@ -60,7 +60,8 @@ SUBROUTINE SEAICE_CALC_LHS( _RL uIceLHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vIceLHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) -#ifdef SEAICE_ALLOW_JFNK +#if ( defined SEAICE_CGRID \ + && ( defined SEAICE_ALLOW_JFNK || defined SEAICE_ALLOW_KRYLOV ) ) C i,j,bi,bj,k :: loop indices INTEGER i,j,bi,bj INTEGER k @@ -206,7 +207,7 @@ SUBROUTINE SEAICE_CALC_LHS( ENDDO ENDDO -#endif /* SEAICE_ALLOW_JFNK */ +#endif /* SEAICE_CGRID, SEAICE_ALLOW_JFNK and KRYLOV */ RETURN END diff --git a/pkg/seaice/seaice_calc_residual.F b/pkg/seaice/seaice_calc_residual.F index a34846237b..77558f21ed 100644 --- a/pkg/seaice/seaice_calc_residual.F +++ b/pkg/seaice/seaice_calc_residual.F @@ -49,7 +49,7 @@ SUBROUTINE SEAICE_CALC_RESIDUAL( _RL uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) -#ifdef SEAICE_ALLOW_JFNK +#if ( defined SEAICE_CGRID && defined SEAICE_ALLOW_JFNK ) C u/vIceLHS :: left hand side of momentum equations _RL uIceLHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vIceLHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) @@ -141,7 +141,7 @@ SUBROUTINE SEAICE_CALC_RESIDUAL( ENDDO ENDDO -#endif /* SEAICE_ALLOW_JFNK */ +#endif /* SEAICE_CGRID and SEAICE_ALLOW_JFNK */ RETURN END diff --git a/pkg/seaice/seaice_calc_rhs.F b/pkg/seaice/seaice_calc_rhs.F index 9a82c5aed8..2c112e2a87 100644 --- a/pkg/seaice/seaice_calc_rhs.F +++ b/pkg/seaice/seaice_calc_rhs.F @@ -47,7 +47,8 @@ SUBROUTINE SEAICE_CALC_RHS( _RL uIceRHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vIceRHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) -#ifdef SEAICE_ALLOW_JFNK +#if ( defined SEAICE_CGRID \ + && ( defined SEAICE_ALLOW_JFNK || defined SEAICE_ALLOW_KRYLOV ) ) C i,j,bi,bj,kSrf :: loop indices INTEGER i,j,bi,bj INTEGER kSrf @@ -116,7 +117,7 @@ SUBROUTINE SEAICE_CALC_RHS( ENDDO ENDDO -#endif /* SEAICE_ALLOW_JFNK */ +#endif /* SEAICE_CGRID, SEAICE_ALLOW_JFNK and KRYLOV */ RETURN END diff --git a/pkg/seaice/seaice_calc_strainrates.F b/pkg/seaice/seaice_calc_strainrates.F index 926b6de249..9dfeba61eb 100644 --- a/pkg/seaice/seaice_calc_strainrates.F +++ b/pkg/seaice/seaice_calc_strainrates.F @@ -57,7 +57,6 @@ SUBROUTINE SEAICE_CALC_STRAINRATES( CEOP #ifdef SEAICE_CGRID -#ifdef SEAICE_ALLOW_DYNAMICS C !LOCAL VARIABLES: C === Local variables === C i,j,bi,bj :: Loop counters @@ -199,7 +198,6 @@ SUBROUTINE SEAICE_CALC_STRAINRATES( #endif #endif /* ALLOW_AUTODIFF */ -#endif /* SEAICE_ALLOW_DYNAMICS */ #endif /* SEAICE_CGRID */ RETURN END diff --git a/pkg/seaice/seaice_calc_viscosities.F b/pkg/seaice/seaice_calc_viscosities.F index 467dfd0237..bbacdc1223 100644 --- a/pkg/seaice/seaice_calc_viscosities.F +++ b/pkg/seaice/seaice_calc_viscosities.F @@ -59,7 +59,7 @@ SUBROUTINE SEAICE_CALC_VISCOSITIES( _RL zetaZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) CEndOfInterface -#if ( defined (SEAICE_CGRID) && defined (SEAICE_ALLOW_DYNAMICS) ) +#ifdef SEAICE_CGRID C === Local variables === C i,j,bi,bj - Loop counters C e11, e12, e22 - components of strain rate tensor @@ -477,6 +477,6 @@ SUBROUTINE SEAICE_CALC_VISCOSITIES( ENDDO ENDDO -#endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID */ +#endif /* SEAICE_CGRID */ RETURN END diff --git a/pkg/seaice/seaice_check.F b/pkg/seaice/seaice_check.F index 621e194a0f..ce8fa20de6 100644 --- a/pkg/seaice/seaice_check.F +++ b/pkg/seaice/seaice_check.F @@ -540,11 +540,19 @@ SUBROUTINE SEAICE_CHECK( myThid ) ENDIF #endif -C-- SEAICE_ALLOW_DYNAMICS and SEAICEuseDYNAMICS -#ifndef SEAICE_ALLOW_DYNAMICS +C-- SEAICE_CGRID, SEAICE_BGRID_DYNAMICS and SEAICEuseDYNAMICS +#if ( defined SEAICE_CGRID && defined SEAICE_BGRID_DYNAMICS ) + WRITE(msgBuf,'(A,A)') + & 'SEAICE_CHECK: SEAICE_CGRID and SEAICE_BGRID_DYNAMICS ', + & 'cannot be defined at the same time.' + CALL PRINT_ERROR( msgBuf, myThid ) + errCount = errCount + 1 +#endif +#if !( defined SEAICE_CGRID || defined SEAICE_BGRID_DYNAMICS ) IF (SEAICEuseDYNAMICS) THEN - WRITE(msgBuf,'(A)') - & 'SEAICE_ALLOW_DYNAMICS needed for SEAICEuseDYNAMICS' + WRITE(msgBuf,'(A,A)') + & 'SEAICE_CHECK: either SEAICE_CGRID ', + & 'or SEAICE_BGRID_DYNAMICS needed for SEAICEuseDYNAMICS = T' CALL PRINT_ERROR( msgBuf, myThid ) errCount = errCount + 1 ENDIF @@ -783,7 +791,7 @@ SUBROUTINE SEAICE_CHECK( myThid ) CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid ) ENDIF -#ifdef SEAICE_ALLOW_DYNAMICS +#ifdef SEAICE_CGRID IF ( SEAICEuseDynamics ) THEN C-- Check Overlap size: IF ( SEAICEuseJFNK ) THEN @@ -873,7 +881,7 @@ SUBROUTINE SEAICE_CHECK( myThid ) ENDIF # endif /* SEAICE_ALLOW_EVP */ -#ifdef SEAICE_ALLOW_EVP +# ifdef SEAICE_ALLOW_EVP IF ( SEAICEuseEVP .AND. (SEAICE_eccfr.NE.SEAICE_eccen) ) THEN WRITE(msgBuf,'(2A)') 'SEAICE_CHECK: SEAICEuseEVP = .TRUE., ', & 'so EVP is turned on by setting appropriate' @@ -904,11 +912,11 @@ SUBROUTINE SEAICE_CHECK( myThid ) CALL PRINT_ERROR( msgBuf, myThid ) errCount = errCount + 1 ENDIF -#endif /* SEAICE_ALLOW_EVP */ +# endif /* SEAICE_ALLOW_EVP */ C SEAICEuseDynamics ENDIF -#endif /* SEAICE_ALLOW_DYNAMICS */ +#endif /* SEAICE_CGRID */ #ifndef SEAICE_GLOBAL_3DIAG_SOLVER IF ( SEAICEuseMultiTileSolver ) THEN @@ -1003,8 +1011,7 @@ SUBROUTINE SEAICE_CHECK( myThid ) errCount = errCount + 1 ENDIF -#ifndef SEAICE_CGRID - +#ifdef SEAICE_BGRID_DYNAMICS IF ( i .GE. 1 ) THEN WRITE(msgBuf,'(2A)') 'SEAICE_CHECK: non-default rheologies ', & 'require that SEAICE_CGRID is defined.' @@ -1036,7 +1043,7 @@ SUBROUTINE SEAICE_CHECK( myThid ) CALL PRINT_ERROR( msgBuf, myThid ) C errCount = errCount + 1 ENDIF -#endif /* ndef SEAICE_CGRID */ +#endif /* SEAICE_BGRID_DYNAMICS */ C-- SEAICE_ALLOW_FREEDRIFT and SEAICEuseFREEDRIFT #ifndef SEAICE_ALLOW_FREEDRIFT diff --git a/pkg/seaice/seaice_dynsolver.F b/pkg/seaice/seaice_dynsolver.F index 398301d1a3..851bcdeffe 100644 --- a/pkg/seaice/seaice_dynsolver.F +++ b/pkg/seaice/seaice_dynsolver.F @@ -8,20 +8,19 @@ C !INTERFACE: SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) -C *==========================================================* -C | SUBROUTINE SEAICE_DYNSOLVER -C | o Ice dynamics using LSR solver -C | Zhang and Hibler, JGR, 102, 8691-8702, 1997 -C | or EVP explicit solver by Hunke and Dukowicz, JPO 27, -C | 1849-1867 (1997) -C *==========================================================* -C | written by Martin Losch, March 2006 -C *==========================================================* -C \ev +C !DESCRIPTION: +C *=============================================================* +C | SUBROUTINE SEAICE_DYNSOLVER | +C | C-grid version of ice dynamics using either | +C | o free drift | +C | o LSR solver, Zhang and Hibler, JGR, 102, 8691-8702, 1997 | +C | o Krylov solver, after Lemieux and Tremblay, JGR, 114, 2009 | +C | o JFNK solver, Losch et al., JCP, 257 901-911, 2014 | +C | o EVP solver, Hunke and Dukowicz, JPO 27, 1849-1867 1997 | +C *=============================================================* C !USES: IMPLICIT NONE - C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" @@ -38,7 +37,7 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) # include "tamc.h" #endif -C !INPUT/OUTPUT PARAMETERS: +C !INPUT PARAMETERS: C === Routine arguments === C myTime :: Simulation time C myIter :: Simulation timestep number @@ -46,23 +45,21 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) _RL myTime INTEGER myIter INTEGER myThid -CEOP -#ifdef SEAICE_CGRID C !FUNCTIONS: LOGICAL DIFFERENT_MULTIPLE EXTERNAL DIFFERENT_MULTIPLE -#ifdef ALLOW_DIAGNOSTICS +#if ( defined ALLOW_DIAGNOSTICS && defined SEAICE_CGRID ) LOGICAL DIAGNOSTICS_IS_ON EXTERNAL DIAGNOSTICS_IS_ON -#endif /* ALLOW_DIAGNOSTICS */ +#endif C !LOCAL VARIABLES: C === Local variables === C i,j :: Loop counters C bi,bj :: tile counters INTEGER i, j, bi, bj -#ifdef SEAICE_ALLOW_DYNAMICS +#ifdef SEAICE_CGRID # ifndef ALLOW_AUTODIFF _RL mask_uice, mask_vice # endif @@ -73,7 +70,7 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) C TAUY :: meridional wind stress over seaice at V point _RL TAUX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL TAUY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) -#ifdef ALLOW_DIAGNOSTICS +#if ( defined ALLOW_DIAGNOSTICS && defined SEAICE_CGRID ) # ifdef ALLOW_AUTODIFF _RL strDivX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL strDivY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) @@ -88,13 +85,14 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) LOGICAL diag_SIsigma_isOn, diag_SIshear_isOn LOGICAL diag_SIenpi_isOn, diag_SIenpot_isOn LOGICAL diag_SIpRfric_isOn, diag_SIpSfric_isOn -#endif /* ALLOW_DIAGNOSTICS */ +#endif +CEOP DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx -#ifdef ALLOW_AUTODIFF +#if ( defined ALLOW_AUTODIFF && defined SEAICE_CGRID ) C Following re-initialisation breaks some "artificial" AD dependencies C incured by IF (DIFFERENT_MULTIPLE ... statement PRESS0 (i,j,bi,bj) = SEAICE_strength*HEFF(i,j,bi,bj) @@ -106,11 +104,11 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) uice_fd(i,j,bi,bj)= 0. _d 0 vice_fd(i,j,bi,bj)= 0. _d 0 # endif -# ifdef ALLOW_DIAGNOSTICS +#if ( defined ALLOW_DIAGNOSTICS && defined SEAICE_CGRID ) strDivX(i,j,bi,bj)= 0. _d 0 strDivY(i,j,bi,bj)= 0. _d 0 # endif -#endif /* ALLOW_AUTODIFF */ +#endif /* ALLOW_AUTODIFF and SEAICE_CGRID */ C Always initialise these local variables, needed for TAF, but also C because they are not completely filled in S/R seaice_get_dynforcing TAUX(i,j,bi,bj) = 0. _d 0 @@ -119,12 +117,12 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) ENDDO ENDDO ENDDO -# ifdef ALLOW_AUTODIFF_TAMC +#ifdef ALLOW_AUTODIFF_TAMC CADJ STORE uice = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE vice = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE uicenm1 = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE vicenm1 = comlev1, key=ikey_dynamics, kind=isbyte -# endif /* ALLOW_AUTODIFF_TAMC */ +#endif /* ALLOW_AUTODIFF_TAMC */ C-- interface of dynamics with atmopheric forcing fields (wind/stress) C Call this in each time step so that we can use the surface stress C in S/R seaice_ocean_stress @@ -133,7 +131,7 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) O TAUX, TAUY, I myTime, myIter, myThid ) -#ifdef SEAICE_ALLOW_DYNAMICS +#ifdef SEAICE_CGRID IF ( SEAICEuseDYNAMICS .AND. & DIFFERENT_MULTIPLE(SEAICE_deltaTdyn,myTime,SEAICE_deltaTtherm) & ) THEN @@ -173,7 +171,7 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) ENDDO ENDDO -#ifndef ALLOW_AUTODIFF +# ifndef ALLOW_AUTODIFF IF ( SEAICE_maskRHS ) THEN C dynamic masking of areas with no ice, not recommended C and only kept for testing purposes @@ -203,12 +201,12 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) ENDDO CALL EXCH_UV_XY_RL( seaiceMaskU, seaiceMaskV, .FALSE., myThid ) ENDIF -#endif /* ndef ALLOW_AUTODIFF */ +# endif /* ndef ALLOW_AUTODIFF */ C-- NOW SET UP FORCING FIELDS C initialise fields -#if (defined ALLOW_AUTODIFF && defined SEAICE_ALLOW_EVP) +# if (defined ALLOW_AUTODIFF && defined SEAICE_ALLOW_EVP) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy @@ -219,7 +217,7 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) ENDDO ENDDO ENDDO -#endif +# endif DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) @@ -238,7 +236,7 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) ENDDO ENDDO ENDIF -#ifdef ATMOSPHERIC_LOADING +# ifdef ATMOSPHERIC_LOADING C-- add atmospheric loading and Sea-Ice loading as it is done for phi0surf C-- in S/R external_forcing_surf IF ( usingZCoords ) THEN @@ -263,7 +261,7 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) C The true atmospheric P-loading is not yet implemented for P-coord C (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS). ENDIF -#endif /* ATMOSPHERIC_LOADING */ +# endif /* ATMOSPHERIC_LOADING */ C-- basic forcing by wind stress IF ( SEAICEscaleSurfStress ) THEN DO j=1-OLy+1,sNy+OLy @@ -302,7 +300,7 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) ENDDO ENDDO -#ifdef SEAICE_ALLOW_FREEDRIFT +# ifdef SEAICE_ALLOW_FREEDRIFT IF ( SEAICEuseFREEDRIFT .OR. SEAICEuseEVP & .OR. LSR_mixIniGuess.EQ.0 ) THEN CALL SEAICE_FREEDRIFT( myTime, myIter, myThid ) @@ -321,70 +319,72 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) ENDDO ENDDO ENDIF -#endif /* SEAICE_ALLOW_FREEDRIFT */ +# endif /* SEAICE_ALLOW_FREEDRIFT */ -#ifdef ALLOW_OBCS +# ifdef ALLOW_OBCS IF ( useOBCS ) THEN CALL OBCS_APPLY_UVICE( uIce, vIce, myThid ) ENDIF -#endif /* ALLOW_OBCS */ +# endif -#ifdef SEAICE_ALLOW_EVP -# ifdef ALLOW_AUTODIFF_TAMC +# ifdef SEAICE_ALLOW_EVP +# ifdef ALLOW_AUTODIFF_TAMC CADJ STORE uice = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE vice = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE uicenm1 = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE vicenm1 = comlev1, key=ikey_dynamics, kind=isbyte -# endif /* ALLOW_AUTODIFF_TAMC */ +# endif /* ALLOW_AUTODIFF_TAMC */ IF ( SEAICEuseEVP ) THEN C Elastic-Viscous-Plastic solver, following Hunke (2001) CALL SEAICE_EVP( myTime, myIter, myThid ) ENDIF -#endif /* SEAICE_ALLOW_EVP */ +# endif /* SEAICE_ALLOW_EVP */ IF ( SEAICEuseLSR ) THEN C Picard solver with LSR scheme (Zhang-J/Hibler 1997), ported to a C-grid CALL SEAICE_LSR( myTime, myIter, myThid ) ENDIF -#ifdef SEAICE_ALLOW_KRYLOV -# ifdef ALLOW_AUTODIFF +# ifdef SEAICE_ALLOW_KRYLOV +# ifdef ALLOW_AUTODIFF STOP 'Adjoint does not work with Picard-Krylov solver.' -# else +# else IF ( SEAICEuseKrylov ) THEN C Picard solver with Matrix-free Krylov solver (Lemieux et al. 2008) CALL SEAICE_KRYLOV( myTime, myIter, myThid ) ENDIF -# endif /* ALLOW_AUTODIFF */ -#endif /* SEAICE_ALLOW_KRYLOV */ +# endif /* ALLOW_AUTODIFF */ +# endif /* SEAICE_ALLOW_KRYLOV */ -#ifdef SEAICE_ALLOW_JFNK -# ifdef ALLOW_AUTODIFF +# ifdef SEAICE_ALLOW_JFNK +# ifdef ALLOW_AUTODIFF STOP 'Adjoint does not work with JFNK solver.' -# else +# else IF ( SEAICEuseJFNK ) THEN C Jacobian-free Newton Krylov solver (Lemieux et al. 2010, 2012) CALL SEAICE_JFNK( myTime, myIter, myThid ) ENDIF -# endif /* ALLOW_AUTODIFF */ -#endif /* SEAICE_ALLOW_JFNK */ +# endif /* ALLOW_AUTODIFF */ +# endif /* SEAICE_ALLOW_JFNK */ C End of IF (SEAICEuseDYNAMICS and DIFFERENT_MULTIPLE ... ENDIF -#endif /* SEAICE_ALLOW_DYNAMICS */ +#endif /* SEAICE_CGRID */ C Update ocean surface stress IF ( SEAICEupdateOceanStress ) THEN #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE uice, vice, DWATN = comlev1, key=ikey_dynamics, kind=isbyte +# ifdef SEAICE_CGRID CADJ STORE stressDivergenceX = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE stressDivergenceY = comlev1, key=ikey_dynamics, kind=isbyte +# endif #endif /* ALLOW_AUTODIFF_TAMC */ CALL SEAICE_OCEAN_STRESS ( I TAUX, TAUY, myTime, myIter, myThid ) ENDIF -#if (defined SEAICE_ALLOW_DYNAMICS && defined SEAICE_ALLOW_CLIPVELS) +#ifdef SEAICE_ALLOW_CLIPVELS IF ( SEAICEuseDYNAMICS .AND. SEAICE_clipVelocities) THEN # ifdef ALLOW_AUTODIFF_TAMC CADJ STORE uice = comlev1, key=ikey_dynamics, kind=isbyte @@ -406,9 +406,9 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) ENDDO ENDDO ENDIF -#endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_ALLOW_CLIPVELS */ +#endif /* SEAICE_ALLOW_CLIPVELS */ -#ifdef ALLOW_DIAGNOSTICS +#if ( defined ALLOW_DIAGNOSTICS && defined SEAICE_CGRID ) C diagnostics related to mechanics/dynamics/momentum equations IF ( useDiagnostics .AND. SEAICEuseDYNAMICS ) THEN CALL DIAGNOSTICS_FILL(zeta ,'SIzeta ',0,1,0,1,1,myThid) @@ -704,10 +704,9 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) C bi/bj-loop ENDDO ENDDO -C useDiagnostics +C useDiagnostics & SEAICEuseDYNAMICS ENDIF -#endif /* ALLOW_DIAGNOSTICS */ -#endif /* SEAICE_CGRID */ +#endif /* ALLOW_DIAGNOSTICS and SEAICE_CGRID */ RETURN END diff --git a/pkg/seaice/seaice_evp.F b/pkg/seaice/seaice_evp.F index 20a0c10805..08e2a890f7 100644 --- a/pkg/seaice/seaice_evp.F +++ b/pkg/seaice/seaice_evp.F @@ -51,9 +51,7 @@ SUBROUTINE SEAICE_EVP( myTime, myIter, myThid ) INTEGER myIter INTEGER myThid -#if ( defined (SEAICE_CGRID) && \ - defined (SEAICE_ALLOW_EVP) && \ - defined (SEAICE_ALLOW_DYNAMICS) ) +#if ( defined SEAICE_CGRID && defined SEAICE_ALLOW_EVP ) C === Local variables === C i,j,bi,bj :: Loop counters @@ -962,7 +960,7 @@ SUBROUTINE SEAICE_EVP( myTime, myIter, myThid ) C end of the main time loop ENDDO -#endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID and SEAICE_ALLOW_EVP */ +#endif /* SEAICE_CGRID and SEAICE_ALLOW_EVP */ RETURN END diff --git a/pkg/seaice/seaice_freedrift.F b/pkg/seaice/seaice_freedrift.F index 099df702d9..328058d73e 100644 --- a/pkg/seaice/seaice_freedrift.F +++ b/pkg/seaice/seaice_freedrift.F @@ -27,8 +27,7 @@ SUBROUTINE SEAICE_FREEDRIFT( myTime, myIter, myThid ) INTEGER myThid CEOP -#ifdef SEAICE_ALLOW_FREEDRIFT -#ifdef SEAICE_CGRID +#if ( defined SEAICE_CGRID && defined SEAICE_ALLOW_FREEDRIFT ) C === Local variables === INTEGER i, j, kSrf, bi, bj @@ -172,7 +171,6 @@ SUBROUTINE SEAICE_FREEDRIFT( myTime, myIter, myThid ) ENDDO ENDDO -#endif /* SEAICE_CGRID */ -#endif /* SEAICE_ALLOW_FREEDRIFT */ +#endif /* SEAICE_CGRID and SEAICE_ALLOW_FREEDRIFT */ RETURN END diff --git a/pkg/seaice/seaice_get_dynforcing.F b/pkg/seaice/seaice_get_dynforcing.F index 0c149e20ec..441a4bc571 100644 --- a/pkg/seaice/seaice_get_dynforcing.F +++ b/pkg/seaice/seaice_get_dynforcing.F @@ -67,7 +67,6 @@ SUBROUTINE SEAICE_GET_DYNFORCING( INTEGER myThid CEOP -#ifdef SEAICE_CGRID C !LOCAL VARIABLES: C i,j,bi,bj :: Loop counters C ks :: vertical index of surface layer @@ -278,7 +277,5 @@ SUBROUTINE SEAICE_GET_DYNFORCING( ENDIF #endif /* ALLOW_DIAGNOSTICS */ -#endif /* SEAICE_CGRID */ - RETURN END diff --git a/pkg/seaice/seaice_init_fixed.F b/pkg/seaice/seaice_init_fixed.F index 5f050caef9..246f930c9a 100644 --- a/pkg/seaice/seaice_init_fixed.F +++ b/pkg/seaice/seaice_init_fixed.F @@ -47,7 +47,7 @@ SUBROUTINE SEAICE_INIT_FIXED( myThid ) #endif /* SHORTWAVE_HEATING */ C local copy of surface layer thickness in meters _RL dzSurf -#ifndef SEAICE_CGRID +#ifdef SEAICE_BGRID_DYNAMICS _RL mask_uice #endif @@ -268,7 +268,7 @@ SUBROUTINE SEAICE_INIT_FIXED( myThid ) SIMaskV(i,j,bi,bj) = maskS(i,j,kSrf,bi,bj) ENDDO ENDDO -#ifndef SEAICE_CGRID +#ifdef SEAICE_BGRID_DYNAMICS DO j=1-OLy+1,sNy+OLy DO i=1-OLx+1,sNx+OLx UVM(i,j,bi,bj)=0. _d 0 @@ -277,7 +277,7 @@ SUBROUTINE SEAICE_INIT_FIXED( myThid ) IF(mask_uice.GT.3.5 _d 0) UVM(i,j,bi,bj)=1. _d 0 ENDDO ENDDO -#endif /* ndef SEAICE_CGRID */ +#endif C coefficients for metric terms #ifdef SEAICE_CGRID @@ -334,7 +334,9 @@ SUBROUTINE SEAICE_INIT_FIXED( myThid ) ENDDO ENDDO ENDIF -#else /* not SEAICE_CGRID */ +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx k1AtC(i,j,bi,bj) = 0.0 _d 0 @@ -405,7 +407,30 @@ SUBROUTINE SEAICE_INIT_FIXED( myThid ) ENDDO ENDDO ENDIF -#endif /* not SEAICE_CGRID */ +C-- Choose a proxy level for geostrophic velocity, + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + KGEO(i,j,bi,bj) = 0 + ENDDO + ENDDO + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx +#ifdef SEAICE_BICE_STRESS + KGEO(i,j,bi,bj) = 1 +#else /* SEAICE_BICE_STRESS */ + IF (klowc(i,j,bi,bj) .LT. 2) THEN + KGEO(i,j,bi,bj) = 1 + ELSE + KGEO(i,j,bi,bj) = 2 + DO WHILE ( abs(rC(KGEO(i,j,bi,bj))) .LT. 50.0 _d 0 .AND. + & KGEO(i,j,bi,bj) .LT. (klowc(i,j,bi,bj)-1) ) + KGEO(i,j,bi,bj) = KGEO(i,j,bi,bj) + 1 + ENDDO + ENDIF +#endif /* SEAICE_BICE_STRESS */ + ENDDO + ENDDO +#endif /* SEAICE_BGRID_DYNAMICS */ C- end bi,bj loops ENDDO ENDDO @@ -434,35 +459,5 @@ SUBROUTINE SEAICE_INIT_FIXED( myThid ) ENDIF #endif /* ALLOW_SHELFICE */ -#ifndef SEAICE_CGRID -C-- Choose a proxy level for geostrophic velocity, - DO bj=myByLo(myThid),myByHi(myThid) - DO bi=myBxLo(myThid),myBxHi(myThid) - DO j=1-OLy,sNy+OLy - DO i=1-OLx,sNx+OLx - KGEO(i,j,bi,bj) = 0 - ENDDO - ENDDO - DO j=1-OLy,sNy+OLy - DO i=1-OLx,sNx+OLx -#ifdef SEAICE_BICE_STRESS - KGEO(i,j,bi,bj) = 1 -#else /* SEAICE_BICE_STRESS */ - IF (klowc(i,j,bi,bj) .LT. 2) THEN - KGEO(i,j,bi,bj) = 1 - ELSE - KGEO(i,j,bi,bj) = 2 - DO WHILE ( abs(rC(KGEO(i,j,bi,bj))) .LT. 50.0 _d 0 .AND. - & KGEO(i,j,bi,bj) .LT. (klowc(i,j,bi,bj)-1) ) - KGEO(i,j,bi,bj) = KGEO(i,j,bi,bj) + 1 - ENDDO - ENDIF -#endif /* SEAICE_BICE_STRESS */ - ENDDO - ENDDO - ENDDO - ENDDO -#endif /* ndef SEAICE_CGRID */ - RETURN END diff --git a/pkg/seaice/seaice_init_varia.F b/pkg/seaice/seaice_init_varia.F index 338a09fe90..692808de92 100644 --- a/pkg/seaice/seaice_init_varia.F +++ b/pkg/seaice/seaice_init_varia.F @@ -38,7 +38,6 @@ SUBROUTINE SEAICE_INIT_VARIA( myThid ) INTEGER i, j, bi, bj INTEGER kSrf - _RS mask_uice INTEGER k #ifdef ALLOW_SITRACER INTEGER iTr, jTh @@ -47,6 +46,7 @@ SUBROUTINE SEAICE_INIT_VARIA( myThid ) INTEGER I_obc, J_obc #endif /* ALLOW_OBCS */ #ifdef SEAICE_CGRID + _RS mask_uice _RL recip_tensilDepth #endif @@ -63,72 +63,69 @@ SUBROUTINE SEAICE_INIT_VARIA( myThid ) DO i=1-OLx,sNx+OLx HEFF(i,j,bi,bj)=0. _d 0 AREA(i,j,bi,bj)=0. _d 0 -CToM<<< + HSNOW(i,j,bi,bj) = 0. _d 0 #ifdef SEAICE_ITD DO k=1,nITD AREAITD(i,j,k,bi,bj) =0. _d 0 HEFFITD(i,j,k,bi,bj) =0. _d 0 + HSNOWITD(i,j,k,bi,bj)=0. _d 0 ENDDO #endif -C>>>ToM UICE(i,j,bi,bj)=0. _d 0 VICE(i,j,bi,bj)=0. _d 0 -#ifdef SEAICE_ALLOW_FREEDRIFT - uice_fd(i,j,bi,bj)=0. _d 0 - vice_fd(i,j,bi,bj)=0. _d 0 -#endif C uIceNm1(i,j,bi,bj) = 0. _d 0 vIceNm1(i,j,bi,bj) = 0. _d 0 + DWATN (i,j,bi,bj) = 0. _d 0 +#ifdef SEAICE_CGRID + stressDivergenceX(i,j,bi,bj) = 0. _d 0 + stressDivergenceY(i,j,bi,bj) = 0. _d 0 +# ifdef SEAICE_ALLOW_EVP + seaice_sigma1 (i,j,bi,bj) = 0. _d 0 + seaice_sigma2 (i,j,bi,bj) = 0. _d 0 + seaice_sigma12(i,j,bi,bj) = 0. _d 0 +# endif /* SEAICE_ALLOW_EVP */ +#endif +#if ( defined SEAICE_CGRID || defined SEAICE_BGRID_DYNAMICS ) e11 (i,j,bi,bj) = 0. _d 0 e22 (i,j,bi,bj) = 0. _d 0 e12 (i,j,bi,bj) = 0. _d 0 deltaC (i,j,bi,bj) = 0. _d 0 + PRESS (i,j,bi,bj) = 0. _d 0 ETA (i,j,bi,bj) = 0. _d 0 etaZ (i,j,bi,bj) = 0. _d 0 ZETA (i,j,bi,bj) = 0. _d 0 FORCEX (i,j,bi,bj) = 0. _d 0 FORCEY (i,j,bi,bj) = 0. _d 0 + tensileStrFac(i,j,bi,bj) = 0. _d 0 + PRESS0(i,j,bi,bj) = 0. _d 0 + FORCEX0(i,j,bi,bj) = 0. _d 0 + FORCEY0(i,j,bi,bj) = 0. _d 0 + SEAICE_zMax(i,j,bi,bj) = 0. _d 0 + SEAICE_zMin(i,j,bi,bj) = 0. _d 0 +#endif #ifdef SEAICE_CGRID seaiceMassC(i,j,bi,bj)=0. _d 0 seaiceMassU(i,j,bi,bj)=0. _d 0 seaiceMassV(i,j,bi,bj)=0. _d 0 - stressDivergenceX(i,j,bi,bj) = 0. _d 0 - stressDivergenceY(i,j,bi,bj) = 0. _d 0 -# ifdef SEAICE_ALLOW_EVP - seaice_sigma1 (i,j,bi,bj) = 0. _d 0 - seaice_sigma2 (i,j,bi,bj) = 0. _d 0 - seaice_sigma12(i,j,bi,bj) = 0. _d 0 -# endif /* SEAICE_ALLOW_EVP */ -#else /* SEAICE_CGRID */ - uIceC(i,j,bi,bj) = 0. _d 0 - vIceC(i,j,bi,bj) = 0. _d 0 +# ifdef SEAICE_ALLOW_FREEDRIFT + uice_fd(i,j,bi,bj) = 0. _d 0 + vice_fd(i,j,bi,bj) = 0. _d 0 +# endif +# ifdef SEAICE_ALLOW_BOTTOMDRAG + CbotC(i,j,bi,bj) = 0. _d 0 +# endif /* SEAICE_ALLOW_BOTTOMDRAG */ +#endif /* SEAICE_CGRID */ +#ifdef SEAICE_BGRID_DYNAMICS + uIceB(i,j,bi,bj) = 0. _d 0 + vIceB(i,j,bi,bj) = 0. _d 0 AMASS(i,j,bi,bj) = 0. _d 0 DAIRN(i,j,bi,bj) = 0. _d 0 WINDX(i,j,bi,bj) = 0. _d 0 WINDY(i,j,bi,bj) = 0. _d 0 GWATX(i,j,bi,bj) = 0. _d 0 GWATY(i,j,bi,bj) = 0. _d 0 -#endif /* SEAICE_CGRID */ - DWATN(i,j,bi,bj) = 0. _d 0 -#ifdef SEAICE_ALLOW_BOTTOMDRAG - CbotC(i,j,bi,bj) = 0. _d 0 -#endif /* SEAICE_ALLOW_BOTTOMDRAG */ - PRESS0(i,j,bi,bj) = 0. _d 0 - PRESS (i,j,bi,bj) = 0. _d 0 - FORCEX0(i,j,bi,bj) = 0. _d 0 - FORCEY0(i,j,bi,bj) = 0. _d 0 - HSNOW(i,j,bi,bj) = 0. _d 0 - SEAICE_zMax (i,j,bi,bj) = 0. _d 0 - SEAICE_zMin (i,j,bi,bj) = 0. _d 0 - tensileStrFac(i,j,bi,bj) = 0. _d 0 -CToM<<< -#ifdef SEAICE_ITD - DO k=1,nITD - HSNOWITD(i,j,k,bi,bj)=0. _d 0 - ENDDO -#endif -C>>>ToM +#endif /* SEAICE_BGRID_DYNAMICS */ #ifdef SEAICE_VARIABLE_SALINITY HSALT(i,j,bi,bj) = 0. _d 0 #endif @@ -185,9 +182,9 @@ SUBROUTINE SEAICE_INIT_VARIA( myThid ) C-- Initialize (variable) grid info. As long as we allow masking of C-- velocities outside of ice covered areas (in seaice_dynsolver) C-- we need to re-initialize seaiceMaskU/V here for TAF/TAMC -#ifdef SEAICE_CGRID DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) +#ifdef SEAICE_CGRID DO j=1-OLy+1,sNy+OLy DO i=1-OLx+1,sNx+OLx seaiceMaskU(i,j,bi,bj)= 0.0 _d 0 @@ -198,14 +195,7 @@ SUBROUTINE SEAICE_INIT_VARIA( myThid ) IF(mask_uice.GT.1.5 _d 0) seaiceMaskV(i,j,bi,bj)=1.0 _d 0 ENDDO ENDDO - ENDDO - ENDDO -#endif /* SEAICE_CGRID */ - - DO bj=myByLo(myThid),myByHi(myThid) - DO bi=myBxLo(myThid),myBxHi(myThid) -#ifdef OBCS_UVICE_OLD -#ifdef SEAICE_CGRID +# ifdef OBCS_UVICE_OLD IF (useOBCS) THEN C-- If OBCS is turned on, close southern and western boundaries DO i=1-OLx,sNx+OLx @@ -225,20 +215,21 @@ SUBROUTINE SEAICE_INIT_VARIA( myThid ) ENDIF ENDDO ENDIF +# endif /* OBCS_UVICE_OLD */ #endif /* SEAICE_CGRID */ -#endif /* OBCS_UVICE_OLD */ DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx DO k=1,nITD TICES(i,j,k,bi,bj)=273.0 _d 0 ENDDO -#ifndef SEAICE_CGRID - AMASS (i,j,bi,bj)=1000.0 _d 0 -#else +#ifdef SEAICE_CGRID seaiceMassC(i,j,bi,bj)=1000.0 _d 0 seaiceMassU(i,j,bi,bj)=1000.0 _d 0 seaiceMassV(i,j,bi,bj)=1000.0 _d 0 +#endif +#ifdef SEAICE_BGRID_DYNAMICS + AMASS (i,j,bi,bj)=1000.0 _d 0 #endif ENDDO ENDDO @@ -249,7 +240,8 @@ SUBROUTINE SEAICE_INIT_VARIA( myThid ) C-- Update overlap regions #ifdef SEAICE_CGRID CALL EXCH_UV_XY_RL(seaiceMaskU,seaiceMaskV,.FALSE.,myThid) -#else +#endif +#ifdef SEAICE_BGRID_DYNAMICS _EXCH_XY_RS(UVM, myThid) #endif @@ -262,7 +254,8 @@ SUBROUTINE SEAICE_INIT_VARIA( myThid ) & nIter0, myThid ) CALL PLOT_FIELD_XYRL( seaiceMaskV, 'Current seaiceMaskV', & nIter0, myThid ) -#else +#endif +#ifdef SEAICE_BGRID_DYNAMICS CALL PLOT_FIELD_XYRS( UVM , 'Current UVM ' , & nIter0, myThid ) #endif @@ -449,25 +442,10 @@ SUBROUTINE SEAICE_INIT_VARIA( myThid ) c ENDIF #endif /* ALLOW_OBCS */ -#ifdef SEAICE_ALLOW_JFNK -C Computing this metric cannot be done in S/R SEAICE_INIT_FIXED -C where it belongs, because globalArea is only defined later after -C S/R PACKAGES_INIT_FIXED, so we move this computation here. - CALL SEAICE_MAP_RS2VEC( nVec, rAw, rAs, - & scalarProductMetric, .TRUE., myThid ) - DO bj=myByLo(myThid),myByHi(myThid) - DO bi=myBxLo(myThid),myBxHi(myThid) - DO i=1,nVec - scalarProductMetric(i,1,bi,bj) = - & scalarProductMetric(i,1,bi,bj)/globalArea - ENDDO - ENDDO - ENDDO -#endif /* SEAICE_ALLOW_JFNK */ - C--- Complete initialization DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) +#if ( defined SEAICE_CGRID || defined SEAICE_BGRID_DYNAMICS ) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx ZETA(i,j,bi,bj) = HEFF(i,j,bi,bj)*(1.0 _d 11) @@ -479,6 +457,7 @@ SUBROUTINE SEAICE_INIT_VARIA( myThid ) PRESS0(i,j,bi,bj) = PRESS0(i,j,bi,bj)*HEFFM(i,j,bi,bj) ENDDO ENDDO +#endif IF ( useRealFreshWaterFlux .AND. .NOT.useThSIce ) THEN DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx @@ -509,6 +488,22 @@ SUBROUTINE SEAICE_INIT_VARIA( myThid ) ENDDO ENDDO ENDIF +# if ( defined SEAICE_ALLOW_JFNK || defined SEAICE_ALLOW_KRYLOV ) +C Computing this metric cannot be done in S/R SEAICE_INIT_FIXED +C where it belongs, because globalArea is only defined later after +C S/R PACKAGES_INIT_FIXED, so we move this computation here. + CALL SEAICE_MAP_RS2VEC( nVec, rAw, rAs, + & scalarProductMetric, .TRUE., myThid ) + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + DO i=1,nVec + scalarProductMetric(i,1,bi,bj) = + & scalarProductMetric(i,1,bi,bj)/globalArea + ENDDO + ENDDO + ENDDO +# endif /* SEAICE_ALLOW_JFNK or KRYLOV */ + #endif /* SEAICE_CGRID */ #ifdef ALLOW_AUTODIFF diff --git a/pkg/seaice/seaice_jfnk.F b/pkg/seaice/seaice_jfnk.F index bf08f80fc0..a82a90cefd 100644 --- a/pkg/seaice/seaice_jfnk.F +++ b/pkg/seaice/seaice_jfnk.F @@ -46,7 +46,7 @@ SUBROUTINE SEAICE_JFNK( myTime, myIter, myThid ) INTEGER myIter INTEGER myThid -#ifdef SEAICE_ALLOW_JFNK +#if ( defined SEAICE_CGRID && defined SEAICE_ALLOW_JFNK ) C !FUNCTIONS: LOGICAL DIFFERENT_MULTIPLE EXTERNAL DIFFERENT_MULTIPLE @@ -595,7 +595,7 @@ SUBROUTINE SEAICE_JFNK_UPDATE( C This is the new residual JFNKresidual = resLoc -#endif /* SEAICE_ALLOW_JFNK */ +#endif /* SEAICE_CGRID and SEAICE_ALLOW_JFNK */ RETURN END diff --git a/pkg/seaice/seaice_krylov.F b/pkg/seaice/seaice_krylov.F index f32f338dc9..a95b620fb5 100644 --- a/pkg/seaice/seaice_krylov.F +++ b/pkg/seaice/seaice_krylov.F @@ -44,7 +44,7 @@ SUBROUTINE SEAICE_KRYLOV( myTime, myIter, myThid ) INTEGER myIter INTEGER myThid -#ifdef SEAICE_ALLOW_KRYLOV +#if ( defined SEAICE_CGRID && defined SEAICE_ALLOW_KRYLOV ) C !FUNCTIONS: LOGICAL DIFFERENT_MULTIPLE EXTERNAL DIFFERENT_MULTIPLE @@ -520,7 +520,7 @@ SUBROUTINE SEAICE_KRYLOV( myTime, myIter, myThid ) _END_MASTER( myThid ) ENDIF -#endif /* SEAICE_ALLOW_KRYLOV */ +#endif /* SEAICE_CGRID and SEAICE_ALLOW_KRYLOV */ RETURN END diff --git a/pkg/seaice/seaice_lsr.F b/pkg/seaice/seaice_lsr.F index bef85641f2..34b90fb582 100644 --- a/pkg/seaice/seaice_lsr.F +++ b/pkg/seaice/seaice_lsr.F @@ -66,7 +66,6 @@ SUBROUTINE SEAICE_LSR( myTime, myIter, myThid ) CEOP #ifdef SEAICE_CGRID -#ifdef SEAICE_ALLOW_DYNAMICS C !FUNCTIONS: LOGICAL DIFFERENT_MULTIPLE EXTERNAL DIFFERENT_MULTIPLE @@ -988,9 +987,6 @@ SUBROUTINE SEAICE_LSR( myTime, myIter, myThid ) C endif useHB87StressCoupling ENDIF -#endif /* SEAICE_ALLOW_DYNAMICS */ -#endif /* SEAICE_CGRID */ - RETURN END @@ -1054,9 +1050,6 @@ SUBROUTINE SEAICE_RESIDUAL( INTEGER myThid CEOP -#ifdef SEAICE_CGRID -#ifdef SEAICE_ALLOW_DYNAMICS - C !LOCAL VARIABLES: C === Local variables === C i,j,bi,bj :: Loop counters @@ -1116,9 +1109,6 @@ SUBROUTINE SEAICE_RESIDUAL( C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| -#endif /* SEAICE_ALLOW_DYNAMICS */ -#endif /* SEAICE_CGRID */ - RETURN END @@ -1184,8 +1174,6 @@ SUBROUTINE SEAICE_LSR_CALC_COEFFS( _RL vRt2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) CEOP -#ifdef SEAICE_CGRID -#ifdef SEAICE_ALLOW_DYNAMICS C !LOCAL VARIABLES: C === Local variables === C i,j,bi,bj :: Loop counters @@ -1997,7 +1985,6 @@ SUBROUTINE SEAICE_LSR_TRIDIAGV( ENDDO #endif /* SEAICE_LSR_ZEBRA */ -#endif /* SEAICE_ALLOW_DYNAMICS */ #endif /* SEAICE_CGRID */ RETURN diff --git a/pkg/seaice/seaice_model.F b/pkg/seaice/seaice_model.F index 4b205e90a1..4c4485eb7d 100644 --- a/pkg/seaice/seaice_model.F +++ b/pkg/seaice/seaice_model.F @@ -23,7 +23,8 @@ SUBROUTINE SEAICE_MODEL( myTime, myIter, myThid ) 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 | o The code has been rewritten substantially to use the | +C | MITgcm C-grid, see Losch et al. OM, 33, 129- 144, 2010 | C *===========================================================* IMPLICIT NONE C \ev @@ -64,7 +65,7 @@ SUBROUTINE SEAICE_MODEL( myTime, myIter, myThid ) #ifdef ALLOW_SITRACER INTEGER iTr #endif -#ifndef SEAICE_CGRID +#ifdef SEAICE_BGRID_DYNAMICS _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 @@ -147,7 +148,7 @@ SUBROUTINE SEAICE_MODEL( myTime, myIter, myThid ) CADJ STORE heff = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE hsnow = comlev1, key=ikey_dynamics, kind=isbyte -# if (defined SEAICE_ALLOW_DYNAMICS) && (defined SEAICE_CGRID) +# ifdef SEAICE_CGRID CADJ STORE fu = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE fv = comlev1, key=ikey_dynamics, kind=isbyte # ifdef SEAICE_ALLOW_EVP @@ -155,26 +156,28 @@ SUBROUTINE SEAICE_MODEL( myTime, myIter, myThid ) CADJ STORE seaice_sigma2 = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE seaice_sigma12 = comlev1, key=ikey_dynamics, kind=isbyte # endif /* SEAICE_ALLOW_EVP */ -# endif /* SEAICE_ALLOW_DYNAMICS) && SEAICE_CGRID */ +# endif /* SEAICE_CGRID */ # 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 +C-- Solve ice momentum equations and calculate ocean surface stress. +C The surface stress always needs to be updated, even if neither B- +C or C-grid dynamics are compiled, and SEAICEuseDYNAMICS = F. #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 +#ifdef SEAICE_BGRID_DYNAMICS CALL TIMER_START('DYNSOLVER [SEAICE_MODEL]',myThid) CALL DYNSOLVER ( myTime, myIter, myThid ) CALL TIMER_STOP ('DYNSOLVER [SEAICE_MODEL]',myThid) -#endif /* SEAICE_CGRID */ +#else /* use default C-grid solver */ + CALL TIMER_START('SEAICE_DYNSOLVER [SEAICE_MODEL]',myThid) + CALL SEAICE_DYNSOLVER ( myTime, myIter, myThid ) + CALL TIMER_STOP ('SEAICE_DYNSOLVER [SEAICE_MODEL]',myThid) +#endif C-- Apply ice velocity open boundary conditions #ifdef ALLOW_OBCS @@ -183,8 +186,7 @@ SUBROUTINE SEAICE_MODEL( myTime, myIter, myThid ) # endif /* DISABLE_SEAICE_OBCS */ #endif /* ALLOW_OBCS */ -#if (defined ALLOW_AUTODIFF_TAMC) && (defined SEAICE_ALLOW_DYNAMICS) \ - && (defined SEAICE_CGRID) +#if ( defined ALLOW_AUTODIFF_TAMC && defined SEAICE_CGRID ) CADJ STORE uice = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE vice = comlev1, key=ikey_dynamics, kind=isbyte C Note: Storing u/vice **after** seaice_dynsolver (and obcs_adjust_uvice) @@ -193,7 +195,7 @@ SUBROUTINE SEAICE_MODEL( myTime, myIter, myThid ) C complicated code in the backward integration (see pkg/autodiff), C the extra call would update the ice velocities with the wrong C set of flag values. -#endif /* ALLOW_AUTODIFF_TAMC && SEAICE_ALLOW_DYNAMICS && SEAICE_CGRID */ +#endif #ifdef ALLOW_THSICE IF ( useThSice ) THEN @@ -213,11 +215,13 @@ SUBROUTINE SEAICE_MODEL( myTime, myIter, myThid ) #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 */ +C-- There always needs to be advection, even if neither B- or C-grid +C dynamics are compiled. +#ifdef SEAICE_BGRID_DYNAMICS CALL SEAICE_ADVDIFF( uLoc, vLoc, myTime, myIter, myThid ) -#endif /* SEAICE_CGRID */ +#else /* default C-grid advection */ + CALL SEAICE_ADVDIFF( uIce, vIce, myTime, myIter, myThid ) +#endif ENDIF #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE heff = comlev1, key=ikey_dynamics, kind=isbyte @@ -253,17 +257,17 @@ SUBROUTINE SEAICE_MODEL( myTime, myIter, myThid ) C must call growth after calling advection C because of ugly time level business IF ( usePW79thermodynamics ) THEN -#ifdef SEAICE_USE_GROWTH_ADX -#ifdef ALLOW_DEBUG +# ifdef SEAICE_USE_GROWTH_ADX +# ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL( 'SEAICE_GROWTH_ADX', myThid ) -#endif +# endif CALL SEAICE_GROWTH_ADX( myTime, myIter, myThid ) -#else -#ifdef ALLOW_DEBUG +# else / * SEAICE_USE_GROWTH_ADX */ +# ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL( 'SEAICE_GROWTH', myThid ) -#endif +# endif CALL SEAICE_GROWTH( myTime, myIter, myThid ) -#endif +# endif / * SEAICE_USE_GROWTH_ADX */ ELSE #endif /* DISABLE_SEAICE_GROWTH */ DO bj=myByLo(myThid),myByHi(myThid) diff --git a/pkg/seaice/seaice_monitor.F b/pkg/seaice/seaice_monitor.F index 9e18f471a5..089e8d76c2 100644 --- a/pkg/seaice/seaice_monitor.F +++ b/pkg/seaice/seaice_monitor.F @@ -100,14 +100,12 @@ SUBROUTINE SEAICE_MONITOR( #ifdef SEAICE_CGRID CALL MON_WRITESTATS_RL( 1, UICE, '_uice', & maskInW, maskInW, rAw, drF, dummyRL, myThid ) -#else - CALL MON_WRITESTATS_RL( 1, UICE, '_uice', - & UVM, UVM, rAz, drF, dummyRL, myThid ) -#endif -#ifdef SEAICE_CGRID CALL MON_WRITESTATS_RL( 1, VICE, '_vice', & maskInS, maskInS, rAs, drF, dummyRL, myThid ) -#else +#endif +#ifdef SEAICE_BGRID_DYNAMICS + CALL MON_WRITESTATS_RL( 1, UICE, '_uice', + & UVM, UVM, rAz, drF, dummyRL, myThid ) CALL MON_WRITESTATS_RL( 1, VICE, '_vice', & UVM, UVM, rAz, drF, dummyRL, myThid ) #endif diff --git a/pkg/seaice/seaice_monitor_ad.F b/pkg/seaice/seaice_monitor_ad.F index 170f7e1152..818c9acb13 100644 --- a/pkg/seaice/seaice_monitor_ad.F +++ b/pkg/seaice/seaice_monitor_ad.F @@ -109,14 +109,12 @@ SUBROUTINE ADSEAICE_MONITOR( #ifdef SEAICE_CGRID CALL MON_WRITESTATS_RL( 1, ADUICE, '_aduice', & maskInW, maskInW, rAw, drF, dummyRL, myThid ) -#else - CALL MON_WRITESTATS_RL( 1, ADUICE, '_aduice', - & UVM, UVM, rAz, drF, dummyRL, myThid ) -#endif -#ifdef SEAICE_CGRID CALL MON_WRITESTATS_RL( 1, ADVICE, '_advice', & maskInS, maskInS, rAs, drF, dummyRL, myThid ) -#else +#endif +#ifdef SEAICE_BGRID_DYNAMICS + CALL MON_WRITESTATS_RL( 1, ADUICE, '_aduice', + & UVM, UVM, rAz, drF, dummyRL, myThid ) CALL MON_WRITESTATS_RL( 1, ADVICE, '_advice', & UVM, UVM, rAz, drF, dummyRL, myThid ) #endif diff --git a/pkg/seaice/seaice_ocean_stress.F b/pkg/seaice/seaice_ocean_stress.F index cc7abe5ece..2480649fb8 100644 --- a/pkg/seaice/seaice_ocean_stress.F +++ b/pkg/seaice/seaice_ocean_stress.F @@ -40,7 +40,6 @@ SUBROUTINE SEAICE_OCEAN_STRESS( INTEGER myIter INTEGER myThid -#ifdef SEAICE_CGRID C !LOCAL VARIABLES: C i,j,bi,bj :: Loop counters C kSrf :: vertical index of surface layer @@ -77,7 +76,9 @@ SUBROUTINE SEAICE_OCEAN_STRESS( & * SEAICEstressFactor fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj) & + areaW*windTauX(I,J,bi,bj) +#ifdef SEAICE_CGRID & + stressDivergenceX(I,J,bi,bj) * SEAICEstressFactor +#endif ENDDO ENDDO C This loop separation makes the adjoint code vectorize @@ -87,7 +88,9 @@ SUBROUTINE SEAICE_OCEAN_STRESS( & * SEAICEstressFactor fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj) & + areaS*windTauY(I,J,bi,bj) +#ifdef SEAICE_CGRID & + stressDivergenceY(I,J,bi,bj) * SEAICEstressFactor +#endif ENDDO ENDDO ENDDO @@ -138,7 +141,5 @@ SUBROUTINE SEAICE_OCEAN_STRESS( ENDIF CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid) -#endif /* SEAICE_CGRID */ - RETURN END diff --git a/pkg/seaice/seaice_oceandrag_coeffs.F b/pkg/seaice/seaice_oceandrag_coeffs.F index ae35a61ef4..8386d3c7a2 100644 --- a/pkg/seaice/seaice_oceandrag_coeffs.F +++ b/pkg/seaice/seaice_oceandrag_coeffs.F @@ -51,7 +51,7 @@ SUBROUTINE SEAICE_OCEANDRAG_COEFFS( C CwatC :: drag coefficients _RL CwatC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) -#if ( (defined SEAICE_CGRID) && (defined SEAICE_ALLOW_DYNAMICS) ) +#ifdef SEAICE_CGRID C === local variables === C i,j,bi,bj,ksrf :: loop indices INTEGER i,j,bi,bj @@ -110,7 +110,7 @@ SUBROUTINE SEAICE_OCEANDRAG_COEFFS( ENDDO ENDDO -#endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID */ +#endif /* SEAICE_CGRID */ RETURN END diff --git a/pkg/seaice/seaice_preconditioner.F b/pkg/seaice/seaice_preconditioner.F index 0a280771a3..74de68182d 100644 --- a/pkg/seaice/seaice_preconditioner.F +++ b/pkg/seaice/seaice_preconditioner.F @@ -68,7 +68,8 @@ SUBROUTINE SEAICE_PRECONDITIONER( _RL dvIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) CEOP -#if (defined SEAICE_ALLOW_DYNAMICS && defined SEAICE_ALLOW_JFNK) +#if ( defined SEAICE_CGRID && \ + ( defined SEAICE_ALLOW_JFNK || defined SEAICE_ALLOW_KRYLOV ) ) C !FUNCTIONS: LOGICAL DIFFERENT_MULTIPLE EXTERNAL DIFFERENT_MULTIPLE @@ -633,7 +634,7 @@ SUBROUTINE SEAICE_PRECOND_RHSV ( ENDDO ENDIF -#endif /* SEAICE_ALLOW_JFNK */ +#endif /* SEAICE_CGRID, SEAICE_ALLOW_JFNK and KRYLOV */ RETURN END diff --git a/pkg/seaice/seaice_readparms.F b/pkg/seaice/seaice_readparms.F index 4d8c8c15a4..1975e9210d 100644 --- a/pkg/seaice/seaice_readparms.F +++ b/pkg/seaice/seaice_readparms.F @@ -234,7 +234,7 @@ SUBROUTINE SEAICE_READPARMS( myThid ) _BEGIN_MASTER(myThid) C-- set default sea ice parameters -#ifdef SEAICE_ALLOW_DYNAMICS +#if ( defined SEAICE_CGRID || defined SEAICE_BGRID_DYNAMICS ) SEAICEuseDYNAMICS = .TRUE. #else SEAICEuseDYNAMICS = .FALSE. diff --git a/verification/1D_ocean_ice_column/code/SEAICE_OPTIONS.h b/verification/1D_ocean_ice_column/code/SEAICE_OPTIONS.h index 4fa80e74d7..b52657adb4 100644 --- a/verification/1D_ocean_ice_column/code/SEAICE_OPTIONS.h +++ b/verification/1D_ocean_ice_column/code/SEAICE_OPTIONS.h @@ -1,3 +1,8 @@ +#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. @@ -6,24 +11,12 @@ C | Use this file for selecting options within the sea ice C | package. C *==========================================================* -#ifndef SEAICE_OPTIONS_H -#define SEAICE_OPTIONS_H -#include "PACKAGES_CONFIG.h" -#include "CPP_OPTIONS.h" - #ifdef ALLOW_SEAICE -C Package-specific Options & Macros go here +C--- Package-specific Options & Macros go here C-- Write "text-plots" of certain fields in STDOUT for debugging. #undef SEAICE_DEBUG -C-- Allow sea-ice dynamic code. -C This option is provided to allow use of TAMC -C on the thermodynamics component of the code only. -C Sea-ice dynamics can also be turned off at runtime -C using variable SEAICEuseDYNAMICS. -#define SEAICE_ALLOW_DYNAMICS - 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 @@ -42,8 +35,6 @@ C-- run with sea Ice Thickness Distribution (ITD); C set number of categories (nITD) in SEAICE_SIZE.h #undef SEAICE_ITD -#undef SEAICE_MODIFY_GROWTH_ADJ - C-- Since the missing sublimation term is now included C this flag is needed for backward compatibility #define SEAICE_DISABLE_SUBLIM @@ -57,66 +48,116 @@ C seaice with a fraction (=SEAICE_saltFrac) of freezing seawater salinity. C- Note: SItracer also offers an alternative way to handle variable salinity. #define SEAICE_VARIABLE_SALINITY -C-- Tracers of ice and/or ice cover. -#undef ALLOW_SITRACER -#ifdef ALLOW_SITRACER -C-- To try avoid 'spontaneous generation' of tracer maxima by advdiff. -# define ALLOW_SITRACER_ADVCAP -#endif - -C-- Enable grease ice parameterization -C The grease ice parameterization delays formation of solid -C sea ice from frazil ice by a time constant and provides a -C dynamic calculation of the initial solid sea ice thickness -C HO as a function of winds, currents and available grease ice -C volume. Grease ice does not significantly reduce heat loss -C from the ocean in winter and area covered by grease is thus -C handled like open water. -C (For details see Smedsrud and Martin, 2014, Ann.Glac.) +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 -C to yield grease ice volume. Additionally, the actual grease ice -C layer thickness (diagnostic SIgrsLT) can be saved. +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-- grease ice uses SItracer: + +C-- Tracers of ice and/or ice cover. #ifdef SEAICE_GREASE +C SEAICE_GREASE code requires to define ALLOW_SITRACER # define ALLOW_SITRACER -# define ALLOW_SITRACER_ADVCAP +#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-- Allow sea-ice dynamic code. These options are 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). 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 +C discretization should still work but it is not longer actively +C tested and supported. Define this flag to compile it. It cannot be +C defined together with SEAICE_CGRID +#define SEAICE_BGRID_DYNAMICS + +C-- The following flag should always be set in order to use C the +C-- operational C-grid discretization. +#undef SEAICE_CGRID -C-- Only for the C-grid version it is possible to #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 # define 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- 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 @@ -124,54 +165,108 @@ 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-- By default for B-grid dynamics solver wind stress under sea-ice is + +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS +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- 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- 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-- When set limit the Ice-Loading to mass of 1/5 of Surface ocean grid-box +#endif /* SEAICE_BGRID_DYNAMICS */ + +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- 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- 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 +C This options excludes more complex physics such +C as sublimation, ITD, and frazil. +#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- 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) +C- Sea-ice volume (requires pkg/cost) #undef ALLOW_COST_ICE -c >>> Sea-ice misfit to obs (requires pkg/cost and ecco) -#undef ALLOW_SEAICE_COST_SMR_AREA +#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 */ - -CEH3 ;;; Local Variables: *** -CEH3 ;;; mode:fortran *** -CEH3 ;;; End: *** diff --git a/verification/1D_ocean_ice_column/code_ad/SEAICE_OPTIONS.h b/verification/1D_ocean_ice_column/code_ad/SEAICE_OPTIONS.h index c673406748..e73d735ef1 100644 --- a/verification/1D_ocean_ice_column/code_ad/SEAICE_OPTIONS.h +++ b/verification/1D_ocean_ice_column/code_ad/SEAICE_OPTIONS.h @@ -1,3 +1,8 @@ +#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. @@ -6,29 +11,12 @@ C | Use this file for selecting options within the sea ice C | package. C *==========================================================* -#ifndef SEAICE_OPTIONS_H -#define SEAICE_OPTIONS_H -#include "PACKAGES_CONFIG.h" -#include "CPP_OPTIONS.h" - #ifdef ALLOW_SEAICE -C Package-specific Options & Macros go here - -C Moved here from obsolete ECCO_CPPOPTIONS.h -cph >>>>>> !!!!!! SPECIAL SEAICE FLAG FOR TESTING !!!!!! <<<<<< -#define SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING -cph >>>>>> !!!!!! SPECIAL SEAICE FLAG FOR TESTING !!!!!! <<<<<< +C--- Package-specific Options & Macros go here C-- Write "text-plots" of certain fields in STDOUT for debugging. #undef SEAICE_DEBUG -C-- Allow sea-ice dynamic code. -C This option is provided to allow use of TAMC -C on the thermodynamics component of the code only. -C Sea-ice dynamics can also be turned off at runtime -C using variable SEAICEuseDYNAMICS. -#undef SEAICE_ALLOW_DYNAMICS - 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 @@ -47,8 +35,6 @@ C-- run with sea Ice Thickness Distribution (ITD); C set number of categories (nITD) in SEAICE_SIZE.h #undef SEAICE_ITD -#undef SEAICE_MODIFY_GROWTH_ADJ - C-- Since the missing sublimation term is now included C this flag is needed for backward compatibility #define SEAICE_DISABLE_SUBLIM @@ -62,66 +48,116 @@ C seaice with a fraction (=SEAICE_saltFrac) of freezing seawater salinity. C- Note: SItracer also offers an alternative way to handle variable salinity. #define SEAICE_VARIABLE_SALINITY -C-- Tracers of ice and/or ice cover. -#undef ALLOW_SITRACER -#ifdef ALLOW_SITRACER -C-- To try avoid 'spontaneous generation' of tracer maxima by advdiff. -# define ALLOW_SITRACER_ADVCAP -#endif - -C-- Enable grease ice parameterization -C The grease ice parameterization delays formation of solid -C sea ice from frazil ice by a time constant and provides a -C dynamic calculation of the initial solid sea ice thickness -C HO as a function of winds, currents and available grease ice -C volume. Grease ice does not significantly reduce heat loss -C from the ocean in winter and area covered by grease is thus -C handled like open water. -C (For details see Smedsrud and Martin, 2014, Ann.Glac.) +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 -C to yield grease ice volume. Additionally, the actual grease ice -C layer thickness (diagnostic SIgrsLT) can be saved. +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-- grease ice uses SItracer: + +C-- Tracers of ice and/or ice cover. #ifdef SEAICE_GREASE +C SEAICE_GREASE code requires to define ALLOW_SITRACER # define ALLOW_SITRACER -# define ALLOW_SITRACER_ADVCAP +#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-- Allow sea-ice dynamic code. These options are 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). 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 +C discretization should still work but it is not longer actively +C tested and supported. Define this flag to compile it. It cannot be +C defined together with SEAICE_CGRID +#undef SEAICE_BGRID_DYNAMICS + +C-- The following flag should always be set in order to use C the +C-- operational C-grid discretization. +#undef SEAICE_CGRID -C-- Only for the C-grid version it is possible to #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- 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 @@ -129,52 +165,108 @@ 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-- By default for B-grid dynamics solver wind stress under sea-ice is + +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS +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- 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- 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-- When set limit the Ice-Loading to mass of 1/5 of Surface ocean grid-box +#endif /* SEAICE_BGRID_DYNAMICS */ + +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- 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- 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 +#define 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 +C This options excludes more complex physics such +C as sublimation, ITD, and frazil. +#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- 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) +C- Sea-ice volume (requires pkg/cost) #define 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 */ - -CEH3 ;;; Local Variables: *** -CEH3 ;;; mode:fortran *** -CEH3 ;;; End: *** diff --git a/verification/cpl_aim+ocn/code_ocn/SEAICE_OPTIONS.h b/verification/cpl_aim+ocn/code_ocn/SEAICE_OPTIONS.h index 5d6956b933..4cc1e1b511 100644 --- a/verification/cpl_aim+ocn/code_ocn/SEAICE_OPTIONS.h +++ b/verification/cpl_aim+ocn/code_ocn/SEAICE_OPTIONS.h @@ -1,3 +1,8 @@ +#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. @@ -6,24 +11,12 @@ C | Use this file for selecting options within the sea ice C | package. C *==========================================================* -#ifndef SEAICE_OPTIONS_H -#define SEAICE_OPTIONS_H -#include "PACKAGES_CONFIG.h" -#include "CPP_OPTIONS.h" - #ifdef ALLOW_SEAICE -C Package-specific Options & Macros go here +C--- Package-specific Options & Macros go here C-- Write "text-plots" of certain fields in STDOUT for debugging. #undef SEAICE_DEBUG -C-- Allow sea-ice dynamic code. -C This option is provided to allow use of TAMC -C on the thermodynamics component of the code only. -C Sea-ice dynamics can also be turned off at runtime -C using variable SEAICEuseDYNAMICS. -#define SEAICE_ALLOW_DYNAMICS - 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 @@ -55,76 +48,116 @@ 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-- Tracers of ice and/or ice cover. -#undef ALLOW_SITRACER -#ifdef ALLOW_SITRACER -C-- To try avoid 'spontaneous generation' of tracer maxima by advdiff. -# define ALLOW_SITRACER_ADVCAP -#endif - -C-- Enable grease ice parameterization -C The grease ice parameterization delays formation of solid -C sea ice from frazil ice by a time constant and provides a -C dynamic calculation of the initial solid sea ice thickness -C HO as a function of winds, currents and available grease ice -C volume. Grease ice does not significantly reduce heat loss -C from the ocean in winter and area covered by grease is thus -C handled like open water. -C (For details see Smedsrud and Martin, 2014, Ann.Glac.) +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 -C to yield grease ice volume. Additionally, the actual grease ice -C layer thickness (diagnostic SIgrsLT) can be saved. +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-- grease ice uses SItracer: + +C-- Tracers of ice and/or ice cover. #ifdef SEAICE_GREASE +C SEAICE_GREASE code requires to define ALLOW_SITRACER # define ALLOW_SITRACER -# define ALLOW_SITRACER_ADVCAP +#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-- Allow sea-ice dynamic code. These options are 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). 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. +C discretization should still work but it is not longer actively +C tested and supported. Define this flag to compile it. It cannot be +C defined together with SEAICE_CGRID +#undef SEAICE_BGRID_DYNAMICS + +C-- The following flag should always be set in order to use C the +C-- operational C-grid discretization. #define SEAICE_CGRID -C-- Only for the C-grid version it is possible to #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 # define SEAICE_ALLOW_JFNK + C enable Krylov code by defining the following flag # define 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 # undef SEAICE_ALLOW_EVP # ifdef SEAICE_ALLOW_EVP -C-- When set use SEAICE_zetaMin and SEAICE_evpDampC to limit viscosities +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 # define 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 @@ -132,58 +165,108 @@ 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-- By default for B-grid dynamics solver wind stress under sea-ice is + +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS +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- 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- 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-- When set limit the Ice-Loading to mass of 1/5 of Surface ocean grid-box +#endif /* SEAICE_BGRID_DYNAMICS */ + +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- When set use SEAICE_clipVelocties = .true., to clip U/VICE at 40cm/s, C not recommended #define SEAICE_ALLOW_CLIPVELS -C-- When set cap the sublimation latent heat flux in solve4temp according + +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 +C This options excludes more complex physics such +C as sublimation, ITD, and frazil. #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- 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) +C- Sea-ice volume (requires pkg/cost) #undef ALLOW_COST_ICE -c >>> Sea-ice misfit to obs (requires pkg/cost and ecco) -#undef ALLOW_SEAICE_COST_SMR_AREA +#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 */ - -CEH3 ;;; Local Variables: *** -CEH3 ;;; mode:fortran *** -CEH3 ;;; End: *** diff --git a/verification/global_ocean.cs32x15/code/SEAICE_OPTIONS.h b/verification/global_ocean.cs32x15/code/SEAICE_OPTIONS.h index 25d7cc5b31..6e9eed5559 100644 --- a/verification/global_ocean.cs32x15/code/SEAICE_OPTIONS.h +++ b/verification/global_ocean.cs32x15/code/SEAICE_OPTIONS.h @@ -1,3 +1,8 @@ +#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. @@ -6,24 +11,12 @@ C | Use this file for selecting options within the sea ice C | package. C *==========================================================* -#ifndef SEAICE_OPTIONS_H -#define SEAICE_OPTIONS_H -#include "PACKAGES_CONFIG.h" -#include "CPP_OPTIONS.h" - #ifdef ALLOW_SEAICE -C Package-specific Options & Macros go here +C--- Package-specific Options & Macros go here C-- Write "text-plots" of certain fields in STDOUT for debugging. #undef SEAICE_DEBUG -C-- Allow sea-ice dynamic code. -C This option is provided to allow use of TAMC -C on the thermodynamics component of the code only. -C Sea-ice dynamics can also be turned off at runtime -C using variable SEAICEuseDYNAMICS. -#define SEAICE_ALLOW_DYNAMICS - 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 @@ -55,66 +48,116 @@ 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-- Tracers of ice and/or ice cover. -#undef ALLOW_SITRACER -#ifdef ALLOW_SITRACER -C-- To try avoid 'spontaneous generation' of tracer maxima by advdiff. -# define ALLOW_SITRACER_ADVCAP -#endif - -C-- Enable grease ice parameterization -C The grease ice parameterization delays formation of solid -C sea ice from frazil ice by a time constant and provides a -C dynamic calculation of the initial solid sea ice thickness -C HO as a function of winds, currents and available grease ice -C volume. Grease ice does not significantly reduce heat loss -C from the ocean in winter and area covered by grease is thus -C handled like open water. -C (For details see Smedsrud and Martin, 2014, Ann.Glac.) +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 -C to yield grease ice volume. Additionally, the actual grease ice -C layer thickness (diagnostic SIgrsLT) can be saved. +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-- grease ice uses SItracer: + +C-- Tracers of ice and/or ice cover. #ifdef SEAICE_GREASE +C SEAICE_GREASE code requires to define ALLOW_SITRACER # define ALLOW_SITRACER -# define ALLOW_SITRACER_ADVCAP +#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-- Allow sea-ice dynamic code. These options are 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). 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. +C discretization should still work but it is not longer actively +C tested and supported. Define this flag to compile it. It cannot be +C defined together with SEAICE_CGRID +#undef SEAICE_BGRID_DYNAMICS + +C-- The following flag should always be set in order to use C the +C-- operational C-grid discretization. #define SEAICE_CGRID -C-- Only for the C-grid version it is possible to #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 # define SEAICE_ALLOW_JFNK + C enable Krylov code by defining the following flag # define 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 # undef SEAICE_ALLOW_EVP # ifdef SEAICE_ALLOW_EVP -C-- When set use SEAICE_zetaMin and SEAICE_evpDampC to limit viscosities +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 @@ -122,54 +165,108 @@ 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-- By default for B-grid dynamics solver wind stress under sea-ice is + +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS +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- 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- 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-- When set limit the Ice-Loading to mass of 1/5 of Surface ocean grid-box +#endif /* SEAICE_BGRID_DYNAMICS */ + +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- When set use SEAICE_clipVelocties = .true., to clip U/VICE at 40cm/s, C not recommended #define SEAICE_ALLOW_CLIPVELS -C-- When set cap the sublimation latent heat flux in solve4temp according + +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 +C This options excludes more complex physics such +C as sublimation, ITD, and frazil. +#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- 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) +C- Sea-ice volume (requires pkg/cost) #undef ALLOW_COST_ICE -c >>> Sea-ice misfit to obs (requires pkg/cost and ecco) -#undef ALLOW_SEAICE_COST_SMR_AREA +#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 */ - -CEH3 ;;; Local Variables: *** -CEH3 ;;; mode:fortran *** -CEH3 ;;; End: *** diff --git a/verification/global_ocean.cs32x15/code_ad/SEAICE_OPTIONS.h b/verification/global_ocean.cs32x15/code_ad/SEAICE_OPTIONS.h index e6ded4d6e7..8a3b7defc3 100644 --- a/verification/global_ocean.cs32x15/code_ad/SEAICE_OPTIONS.h +++ b/verification/global_ocean.cs32x15/code_ad/SEAICE_OPTIONS.h @@ -1,3 +1,8 @@ +#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. @@ -6,24 +11,12 @@ C | Use this file for selecting options within the sea ice C | package. C *==========================================================* -#ifndef SEAICE_OPTIONS_H -#define SEAICE_OPTIONS_H -#include "PACKAGES_CONFIG.h" -#include "CPP_OPTIONS.h" - #ifdef ALLOW_SEAICE -C Package-specific Options & Macros go here +C--- Package-specific Options & Macros go here C-- Write "text-plots" of certain fields in STDOUT for debugging. #undef SEAICE_DEBUG -C-- Allow sea-ice dynamic code. -C This option is provided to allow use of TAMC -C on the thermodynamics component of the code only. -C Sea-ice dynamics can also be turned off at runtime -C using variable SEAICEuseDYNAMICS. -#define SEAICE_ALLOW_DYNAMICS - 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 @@ -55,66 +48,116 @@ 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-- Tracers of ice and/or ice cover. -#undef ALLOW_SITRACER -#ifdef ALLOW_SITRACER -C-- To try avoid 'spontaneous generation' of tracer maxima by advdiff. -# define ALLOW_SITRACER_ADVCAP -#endif - -C-- Enable grease ice parameterization -C The grease ice parameterization delays formation of solid -C sea ice from frazil ice by a time constant and provides a -C dynamic calculation of the initial solid sea ice thickness -C HO as a function of winds, currents and available grease ice -C volume. Grease ice does not significantly reduce heat loss -C from the ocean in winter and area covered by grease is thus -C handled like open water. -C (For details see Smedsrud and Martin, 2014, Ann.Glac.) +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 -C to yield grease ice volume. Additionally, the actual grease ice -C layer thickness (diagnostic SIgrsLT) can be saved. +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-- grease ice uses SItracer: + +C-- Tracers of ice and/or ice cover. #ifdef SEAICE_GREASE +C SEAICE_GREASE code requires to define ALLOW_SITRACER # define ALLOW_SITRACER -# define ALLOW_SITRACER_ADVCAP +#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-- Allow sea-ice dynamic code. These options are 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). 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. +C discretization should still work but it is not longer actively +C tested and supported. Define this flag to compile it. It cannot be +C defined together with SEAICE_CGRID +#undef SEAICE_BGRID_DYNAMICS + +C-- The following flag should always be set in order to use C the +C-- operational C-grid discretization. #define SEAICE_CGRID -C-- Only for the C-grid version it is possible to #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- 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 @@ -122,52 +165,108 @@ 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 # define SEAICE_ALLOW_BOTTOMDRAG -#else /* not SEAICE_CGRID, but old B-grid */ -C-- By default for B-grid dynamics solver wind stress under sea-ice is + +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS +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- 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- 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-- When set limit the Ice-Loading to mass of 1/5 of Surface ocean grid-box +#endif /* SEAICE_BGRID_DYNAMICS */ + +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 + +C- When set use SEAICE_clipVelocties = .true., to clip U/VICE at 40cm/s, +C not recommended #define SEAICE_ALLOW_CLIPVELS -C-- When set cap the sublimation latent heat flux in solve4temp according + +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 +C This options excludes more complex physics such +C as sublimation, ITD, and frazil. +#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- Do not compile/use seaice-related obcs code when using obcs. +#undef DISABLE_SEAICE_OBCS + C-- Enable free drift code #define SEAICE_ALLOW_FREEDRIFT C-- pkg/seaice cost functions compile flags -c >>> Sea-ice volume (requires pkg/cost) +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 */ - -CEH3 ;;; Local Variables: *** -CEH3 ;;; mode:fortran *** -CEH3 ;;; End: *** diff --git a/verification/global_ocean.cs32x15/code_tap/SEAICE_OPTIONS.h b/verification/global_ocean.cs32x15/code_tap/SEAICE_OPTIONS.h index 88dcc7134d..8a3b7defc3 100644 --- a/verification/global_ocean.cs32x15/code_tap/SEAICE_OPTIONS.h +++ b/verification/global_ocean.cs32x15/code_tap/SEAICE_OPTIONS.h @@ -1,3 +1,8 @@ +#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. @@ -6,24 +11,12 @@ C | Use this file for selecting options within the sea ice C | package. C *==========================================================* -#ifndef SEAICE_OPTIONS_H -#define SEAICE_OPTIONS_H -#include "PACKAGES_CONFIG.h" -#include "CPP_OPTIONS.h" - #ifdef ALLOW_SEAICE -C Package-specific Options & Macros go here +C--- Package-specific Options & Macros go here C-- Write "text-plots" of certain fields in STDOUT for debugging. #undef SEAICE_DEBUG -C-- Allow sea-ice dynamic code. -C This option is provided to allow use of TAMC -C on the thermodynamics component of the code only. -C Sea-ice dynamics can also be turned off at runtime -C using variable SEAICEuseDYNAMICS. -#define SEAICE_ALLOW_DYNAMICS - 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 @@ -55,66 +48,116 @@ 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-- Tracers of ice and/or ice cover. -#undef ALLOW_SITRACER -#ifdef ALLOW_SITRACER -C-- To try avoid 'spontaneous generation' of tracer maxima by advdiff. -# define ALLOW_SITRACER_ADVCAP -#endif - -C-- Enable grease ice parameterization -C The grease ice parameterization delays formation of solid -C sea ice from frazil ice by a time constant and provides a -C dynamic calculation of the initial solid sea ice thickness -C HO as a function of winds, currents and available grease ice -C volume. Grease ice does not significantly reduce heat loss -C from the ocean in winter and area covered by grease is thus -C handled like open water. -C (For details see Smedsrud and Martin, 2014, Ann.Glac.) +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 -C to yield grease ice volume. Additionally, the actual grease ice -C layer thickness (diagnostic SIgrsLT) can be saved. +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-- grease ice uses SItracer: + +C-- Tracers of ice and/or ice cover. #ifdef SEAICE_GREASE +C SEAICE_GREASE code requires to define ALLOW_SITRACER # define ALLOW_SITRACER -# define ALLOW_SITRACER_ADVCAP +#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-- Allow sea-ice dynamic code. These options are 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). 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. +C discretization should still work but it is not longer actively +C tested and supported. Define this flag to compile it. It cannot be +C defined together with SEAICE_CGRID +#undef SEAICE_BGRID_DYNAMICS + +C-- The following flag should always be set in order to use C the +C-- operational C-grid discretization. #define SEAICE_CGRID -C-- Only for the C-grid version it is possible to #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- 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 @@ -122,54 +165,108 @@ 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 # define SEAICE_ALLOW_BOTTOMDRAG -#else /* not SEAICE_CGRID, but old B-grid */ -C-- By default for B-grid dynamics solver wind stress under sea-ice is + +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS +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- 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- 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-- When set limit the Ice-Loading to mass of 1/5 of Surface ocean grid-box +#endif /* SEAICE_BGRID_DYNAMICS */ + +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 + +C- When set use SEAICE_clipVelocties = .true., to clip U/VICE at 40cm/s, +C not recommended #define SEAICE_ALLOW_CLIPVELS -C-- When set cap the sublimation latent heat flux in solve4temp according + +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 +C This options excludes more complex physics such +C as sublimation, ITD, and frazil. +#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- Do not compile/use seaice-related obcs code when using obcs. +#undef DISABLE_SEAICE_OBCS + C-- Enable free drift code #define SEAICE_ALLOW_FREEDRIFT C-- pkg/seaice cost functions compile flags -c >>> Sea-ice volume (requires pkg/cost) +C- Sea-ice volume (requires pkg/cost) #undef ALLOW_COST_ICE -c >>> Sea-ice misfit to obs (requires pkg/cost and ecco) -#undef ALLOW_SEAICE_COST_SMR_AREA +#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 */ - -CEH3 ;;; Local Variables: *** -CEH3 ;;; mode:fortran *** -CEH3 ;;; End: *** diff --git a/verification/lab_sea/code/SEAICE_OPTIONS.h b/verification/lab_sea/code/SEAICE_OPTIONS.h index 1cb62dd9c9..447a1eef6e 100644 --- a/verification/lab_sea/code/SEAICE_OPTIONS.h +++ b/verification/lab_sea/code/SEAICE_OPTIONS.h @@ -1,3 +1,8 @@ +#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. @@ -6,24 +11,12 @@ C | Use this file for selecting options within the sea ice C | package. C *==========================================================* -#ifndef SEAICE_OPTIONS_H -#define SEAICE_OPTIONS_H -#include "PACKAGES_CONFIG.h" -#include "CPP_OPTIONS.h" - #ifdef ALLOW_SEAICE -C Package-specific Options & Macros go here +C--- Package-specific Options & Macros go here C-- Write "text-plots" of certain fields in STDOUT for debugging. #undef SEAICE_DEBUG -C-- Allow sea-ice dynamic code. -C This option is provided to allow use of TAMC -C on the thermodynamics component of the code only. -C Sea-ice dynamics can also be turned off at runtime -C using variable SEAICEuseDYNAMICS. -#define SEAICE_ALLOW_DYNAMICS - 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 @@ -55,66 +48,116 @@ 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-- Tracers of ice and/or ice cover. -#define ALLOW_SITRACER -#ifdef ALLOW_SITRACER -C-- To try avoid 'spontaneous generation' of tracer maxima by advdiff. -# define ALLOW_SITRACER_ADVCAP -#endif - -C-- Enable grease ice parameterization -C The grease ice parameterization delays formation of solid -C sea ice from frazil ice by a time constant and provides a -C dynamic calculation of the initial solid sea ice thickness -C HO as a function of winds, currents and available grease ice -C volume. Grease ice does not significantly reduce heat loss -C from the ocean in winter and area covered by grease is thus -C handled like open water. -C (For details see Smedsrud and Martin, 2014, Ann.Glac.) +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 -C to yield grease ice volume. Additionally, the actual grease ice -C layer thickness (diagnostic SIgrsLT) can be saved. +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-- grease ice uses SItracer: + +C-- Tracers of ice and/or ice cover. #ifdef SEAICE_GREASE +C SEAICE_GREASE code requires to define ALLOW_SITRACER +# define ALLOW_SITRACER +#else # define ALLOW_SITRACER -# define ALLOW_SITRACER_ADVCAP #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-- Allow sea-ice dynamic code. These options are 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). 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. +C discretization should still work but it is not longer actively +C tested and supported. Define this flag to compile it. It cannot be +C defined together with SEAICE_CGRID +#undef SEAICE_BGRID_DYNAMICS + +C-- The following flag should always be set in order to use C the +C-- operational C-grid discretization. #define SEAICE_CGRID -C-- Only for the C-grid version it is possible to #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- 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 @@ -122,54 +165,108 @@ 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. # define 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 # define SEAICE_ALLOW_BOTTOMDRAG -#else /* not SEAICE_CGRID, but old B-grid */ -C-- By default for B-grid dynamics solver wind stress under sea-ice is + +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS +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- 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- 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-- When set limit the Ice-Loading to mass of 1/5 of Surface ocean grid-box +#endif /* SEAICE_BGRID_DYNAMICS */ + +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- 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- 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 +C This options excludes more complex physics such +C as sublimation, ITD, and frazil. +#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- Do not compile/use seaice-related obcs code when using obcs. +#undef DISABLE_SEAICE_OBCS + C-- Enable free drift code #define SEAICE_ALLOW_FREEDRIFT C-- pkg/seaice cost functions compile flags -c >>> Sea-ice volume (requires pkg/cost) +C- Sea-ice volume (requires pkg/cost) #undef ALLOW_COST_ICE -c >>> Sea-ice misfit to obs (requires pkg/cost and ecco) -#undef ALLOW_SEAICE_COST_SMR_AREA +#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 */ - -CEH3 ;;; Local Variables: *** -CEH3 ;;; mode:fortran *** -CEH3 ;;; End: *** diff --git a/verification/lab_sea/code_ad/SEAICE_OPTIONS.h b/verification/lab_sea/code_ad/SEAICE_OPTIONS.h index 6e972776fe..d27c394be0 100644 --- a/verification/lab_sea/code_ad/SEAICE_OPTIONS.h +++ b/verification/lab_sea/code_ad/SEAICE_OPTIONS.h @@ -1,3 +1,8 @@ +#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. @@ -6,24 +11,12 @@ C | Use this file for selecting options within the sea ice C | package. C *==========================================================* -#ifndef SEAICE_OPTIONS_H -#define SEAICE_OPTIONS_H -#include "PACKAGES_CONFIG.h" -#include "CPP_OPTIONS.h" - #ifdef ALLOW_SEAICE -C Package-specific Options & Macros go here +C--- Package-specific Options & Macros go here C-- Write "text-plots" of certain fields in STDOUT for debugging. #undef SEAICE_DEBUG -C-- Allow sea-ice dynamic code. -C This option is provided to allow use of TAMC -C on the thermodynamics component of the code only. -C Sea-ice dynamics can also be turned off at runtime -C using variable SEAICEuseDYNAMICS. -#define SEAICE_ALLOW_DYNAMICS - 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 @@ -55,66 +48,116 @@ C seaice with a fraction (=SEAICE_saltFrac) of freezing seawater salinity. C- Note: SItracer also offers an alternative way to handle variable salinity. #define SEAICE_VARIABLE_SALINITY -C-- Tracers of ice and/or ice cover. -#undef ALLOW_SITRACER -#ifdef ALLOW_SITRACER -C-- To try avoid 'spontaneous generation' of tracer maxima by advdiff. -# define ALLOW_SITRACER_ADVCAP -#endif - -C-- Enable grease ice parameterization -C The grease ice parameterization delays formation of solid -C sea ice from frazil ice by a time constant and provides a -C dynamic calculation of the initial solid sea ice thickness -C HO as a function of winds, currents and available grease ice -C volume. Grease ice does not significantly reduce heat loss -C from the ocean in winter and area covered by grease is thus -C handled like open water. -C (For details see Smedsrud and Martin, 2014, Ann.Glac.) +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 -C to yield grease ice volume. Additionally, the actual grease ice -C layer thickness (diagnostic SIgrsLT) can be saved. +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-- grease ice uses SItracer: + +C-- Tracers of ice and/or ice cover. #ifdef SEAICE_GREASE +C SEAICE_GREASE code requires to define ALLOW_SITRACER # define ALLOW_SITRACER -# define ALLOW_SITRACER_ADVCAP +#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-- Allow sea-ice dynamic code. These options are 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). 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. +C discretization should still work but it is not longer actively +C tested and supported. Define this flag to compile it. It cannot be +C defined together with SEAICE_CGRID +#undef SEAICE_BGRID_DYNAMICS + +C-- The following flag should always be set in order to use C the +C-- operational C-grid discretization. #define SEAICE_CGRID -C-- Only for the C-grid version it is possible to #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- 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 @@ -122,55 +165,108 @@ 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. # define 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. # define 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-- By default for B-grid dynamics solver wind stress under sea-ice is + +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS +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- 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- 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-- When set limit the Ice-Loading to mass of 1/5 of Surface ocean grid-box +#endif /* SEAICE_BGRID_DYNAMICS */ + +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- 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- 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 +C This options excludes more complex physics such +C as sublimation, ITD, and frazil. +#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- Do not compile/use seaice-related obcs code when using obcs. +#undef DISABLE_SEAICE_OBCS + C-- Enable free drift code #define SEAICE_ALLOW_FREEDRIFT C-- pkg/seaice cost functions compile flags -c >>> Sea-ice volume (requires pkg/cost) +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 */ - -CEH3 ;;; Local Variables: *** -CEH3 ;;; mode:fortran *** -CEH3 ;;; End: *** diff --git a/verification/lab_sea/code_tap/SEAICE_OPTIONS.h b/verification/lab_sea/code_tap/SEAICE_OPTIONS.h index 8e1df30deb..d27c394be0 100644 --- a/verification/lab_sea/code_tap/SEAICE_OPTIONS.h +++ b/verification/lab_sea/code_tap/SEAICE_OPTIONS.h @@ -1,3 +1,8 @@ +#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. @@ -6,24 +11,12 @@ C | Use this file for selecting options within the sea ice C | package. C *==========================================================* -#ifndef SEAICE_OPTIONS_H -#define SEAICE_OPTIONS_H -#include "PACKAGES_CONFIG.h" -#include "CPP_OPTIONS.h" - #ifdef ALLOW_SEAICE -C Package-specific Options & Macros go here +C--- Package-specific Options & Macros go here C-- Write "text-plots" of certain fields in STDOUT for debugging. #undef SEAICE_DEBUG -C-- Allow sea-ice dynamic code. -C This option is provided to allow use of TAMC -C on the thermodynamics component of the code only. -C Sea-ice dynamics can also be turned off at runtime -C using variable SEAICEuseDYNAMICS. -#define SEAICE_ALLOW_DYNAMICS - 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 @@ -55,66 +48,116 @@ C seaice with a fraction (=SEAICE_saltFrac) of freezing seawater salinity. C- Note: SItracer also offers an alternative way to handle variable salinity. #define SEAICE_VARIABLE_SALINITY -C-- Tracers of ice and/or ice cover. -#undef ALLOW_SITRACER -#ifdef ALLOW_SITRACER -C-- To try avoid 'spontaneous generation' of tracer maxima by advdiff. -# define ALLOW_SITRACER_ADVCAP -#endif - -C-- Enable grease ice parameterization -C The grease ice parameterization delays formation of solid -C sea ice from frazil ice by a time constant and provides a -C dynamic calculation of the initial solid sea ice thickness -C HO as a function of winds, currents and available grease ice -C volume. Grease ice does not significantly reduce heat loss -C from the ocean in winter and area covered by grease is thus -C handled like open water. -C (For details see Smedsrud and Martin, 2014, Ann.Glac.) +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 -C to yield grease ice volume. Additionally, the actual grease ice -C layer thickness (diagnostic SIgrsLT) can be saved. +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-- grease ice uses SItracer: + +C-- Tracers of ice and/or ice cover. #ifdef SEAICE_GREASE +C SEAICE_GREASE code requires to define ALLOW_SITRACER # define ALLOW_SITRACER -# define ALLOW_SITRACER_ADVCAP +#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-- Allow sea-ice dynamic code. These options are 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). 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. +C discretization should still work but it is not longer actively +C tested and supported. Define this flag to compile it. It cannot be +C defined together with SEAICE_CGRID +#undef SEAICE_BGRID_DYNAMICS + +C-- The following flag should always be set in order to use C the +C-- operational C-grid discretization. #define SEAICE_CGRID -C-- Only for the C-grid version it is possible to #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- 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 @@ -122,62 +165,108 @@ 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. # define SEAICE_LSR_ZEBRA -C This flag is also required for an actual adjoint of seaice_lsr; + +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. # define 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-- By default for B-grid dynamics solver wind stress under sea-ice is + +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS +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- 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- 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-- When set limit the Ice-Loading to mass of 1/5 of Surface ocean grid-box +#endif /* SEAICE_BGRID_DYNAMICS */ + +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- 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- 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 +C This options excludes more complex physics such +C as sublimation, ITD, and frazil. +#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- Do not compile/use seaice-related obcs code when using obcs. +#undef DISABLE_SEAICE_OBCS + C-- Enable free drift code #define SEAICE_ALLOW_FREEDRIFT C-- pkg/seaice cost functions compile flags -c >>> Sea-ice volume (requires pkg/cost) +C- Sea-ice volume (requires pkg/cost) #undef ALLOW_COST_ICE -c >>> Sea-ice misfit to obs (requires pkg/cost and ecco) -c#define ALLOW_SEAICE_COST_SMR_AREA +#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 */ -C-- ssg - This flag turns on the more "correct" seaice themodynamics -C by Arash Bigdeli. -C-- Use the adjointable sea-ice thermodynamic model -C in seaice_growth_adx.F instead of seaice_growth.F -#undef SEAICE_USE_GROWTH_ADX #endif /* ALLOW_SEAICE */ #endif /* SEAICE_OPTIONS_H */ - -CEH3 ;;; Local Variables: *** -CEH3 ;;; mode:fortran *** -CEH3 ;;; End: *** diff --git a/verification/offline_exf_seaice/code/SEAICE_OPTIONS.h b/verification/offline_exf_seaice/code/SEAICE_OPTIONS.h index 852c4ab3cc..8967c1329a 100644 --- a/verification/offline_exf_seaice/code/SEAICE_OPTIONS.h +++ b/verification/offline_exf_seaice/code/SEAICE_OPTIONS.h @@ -1,3 +1,8 @@ +#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. @@ -6,24 +11,12 @@ C | Use this file for selecting options within the sea ice C | package. C *==========================================================* -#ifndef SEAICE_OPTIONS_H -#define SEAICE_OPTIONS_H -#include "PACKAGES_CONFIG.h" -#include "CPP_OPTIONS.h" - #ifdef ALLOW_SEAICE -C Package-specific Options & Macros go here +C--- Package-specific Options & Macros go here C-- Write "text-plots" of certain fields in STDOUT for debugging. #undef SEAICE_DEBUG -C-- Allow sea-ice dynamic code. -C This option is provided to allow use of TAMC -C on the thermodynamics component of the code only. -C Sea-ice dynamics can also be turned off at runtime -C using variable SEAICEuseDYNAMICS. -#define SEAICE_ALLOW_DYNAMICS - 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 @@ -55,79 +48,122 @@ 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-- Tracers of ice and/or ice cover. -#undef ALLOW_SITRACER -#ifdef ALLOW_SITRACER -C-- To try avoid 'spontaneous generation' of tracer maxima by advdiff. -# define ALLOW_SITRACER_ADVCAP -#endif - -C-- Enable grease ice parameterization -C The grease ice parameterization delays formation of solid -C sea ice from frazil ice by a time constant and provides a -C dynamic calculation of the initial solid sea ice thickness -C HO as a function of winds, currents and available grease ice -C volume. Grease ice does not significantly reduce heat loss -C from the ocean in winter and area covered by grease is thus -C handled like open water. -C (For details see Smedsrud and Martin, 2014, Ann.Glac.) +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 -C to yield grease ice volume. Additionally, the actual grease ice -C layer thickness (diagnostic SIgrsLT) can be saved. +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-- grease ice uses SItracer: + +C-- Tracers of ice and/or ice cover. #ifdef SEAICE_GREASE +C SEAICE_GREASE code requires to define ALLOW_SITRACER # define ALLOW_SITRACER -# define ALLOW_SITRACER_ADVCAP +#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-- Allow sea-ice dynamic code. These options are 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). 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. +C discretization should still work but it is not longer actively +C tested and supported. Define this flag to compile it. It cannot be +C defined together with SEAICE_CGRID +#undef SEAICE_BGRID_DYNAMICS + +C-- The following flag should always be set in order to use C the +C-- operational C-grid discretization. #define SEAICE_CGRID -C-- Only for the C-grid version it is possible to #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 # define SEAICE_ALLOW_JFNK + C enable Krylov code by defining the following flag # define SEAICE_ALLOW_KRYLOV -C enable this flag to reproduce old verification results for JFNK + +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 # define SEAICE_PRECOND_EXTRA_EXCHANGE + C enable LSR to use global (multi-tile) tri-diagonal solver # define 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- 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 # define SEAICE_ZETA_SMOOTHREG + +C-- Different yield curves within the VP rheology framework C allow the truncated ellipse rheology (runtime flag SEAICEuseTEM) -C or the modified coulombic rheology # define SEAICE_ALLOW_TEM -C allow the use of the Mohr Coulomb rheology (runtime flag -C SEAICEuseFULLMC) as defined in (Ip 1991) /!\ This is known -C to give unstable results, use with caution -# define SEAICE_ALLOW_FULLMC + +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 +# define SEAICE_ALLOW_MCS + +CMLC or the modified coulombic rheology +CMLC allow the use of the Mohr Coulomb rheology (runtime flag +CMLC SEAICEuseFULLMC) as defined in (Ip 1991) /!\ This is known +CMLC to give unstable results, use with caution +CML# define SEAICE_ALLOW_FULLMC + C allow the use of Mohr Coulomb with elliptical plastic potential C (runtime flag SEAICEuseMCE) # define SEAICE_ALLOW_MCE -C allow the teardrop and parabolic lens rheology (runtime flag -C SEAICEuseTD and SEAICEusePL) + +C allow the teardrop and parabolic lens rheology +C (runtime flag SEAICEuseTD and SEAICEusePL) # define 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 @@ -135,58 +171,108 @@ 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-- By default for B-grid dynamics solver wind stress under sea-ice is + +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS +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- 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- 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-- When set limit the Ice-Loading to mass of 1/5 of Surface ocean grid-box +#endif /* SEAICE_BGRID_DYNAMICS */ + +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- 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- 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 +C This options excludes more complex physics such +C as sublimation, ITD, and frazil. #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- Do not compile/use seaice-related obcs code when using obcs. +#undef DISABLE_SEAICE_OBCS + C-- Enable free drift code #define SEAICE_ALLOW_FREEDRIFT C-- pkg/seaice cost functions compile flags -c >>> Sea-ice volume (requires pkg/cost) +C- Sea-ice volume (requires pkg/cost) #undef ALLOW_COST_ICE -c >>> Sea-ice misfit to obs (requires pkg/cost and ecco) -#undef ALLOW_SEAICE_COST_SMR_AREA +#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 */ - -CEH3 ;;; Local Variables: *** -CEH3 ;;; mode:fortran *** -CEH3 ;;; End: *** diff --git a/verification/offline_exf_seaice/code_ad/SEAICE_OPTIONS.h b/verification/offline_exf_seaice/code_ad/SEAICE_OPTIONS.h index d031682d4e..be18d11b1d 100644 --- a/verification/offline_exf_seaice/code_ad/SEAICE_OPTIONS.h +++ b/verification/offline_exf_seaice/code_ad/SEAICE_OPTIONS.h @@ -1,3 +1,8 @@ +#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. @@ -6,24 +11,12 @@ C | Use this file for selecting options within the sea ice C | package. C *==========================================================* -#ifndef SEAICE_OPTIONS_H -#define SEAICE_OPTIONS_H -#include "PACKAGES_CONFIG.h" -#include "CPP_OPTIONS.h" - #ifdef ALLOW_SEAICE -C Package-specific Options & Macros go here +C--- Package-specific Options & Macros go here C-- Write "text-plots" of certain fields in STDOUT for debugging. #undef SEAICE_DEBUG -C-- Allow sea-ice dynamic code. -C This option is provided to allow use of TAMC -C on the thermodynamics component of the code only. -C Sea-ice dynamics can also be turned off at runtime -C using variable SEAICEuseDYNAMICS. -#define SEAICE_ALLOW_DYNAMICS - 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 @@ -55,66 +48,122 @@ 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-- Tracers of ice and/or ice cover. -#undef ALLOW_SITRACER -#ifdef ALLOW_SITRACER -C-- To try avoid 'spontaneous generation' of tracer maxima by advdiff. -# define ALLOW_SITRACER_ADVCAP -#endif - -C-- Enable grease ice parameterization -C The grease ice parameterization delays formation of solid -C sea ice from frazil ice by a time constant and provides a -C dynamic calculation of the initial solid sea ice thickness -C HO as a function of winds, currents and available grease ice -C volume. Grease ice does not significantly reduce heat loss -C from the ocean in winter and area covered by grease is thus -C handled like open water. -C (For details see Smedsrud and Martin, 2014, Ann.Glac.) +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 -C to yield grease ice volume. Additionally, the actual grease ice -C layer thickness (diagnostic SIgrsLT) can be saved. +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-- grease ice uses SItracer: + +C-- Tracers of ice and/or ice cover. #ifdef SEAICE_GREASE +C SEAICE_GREASE code requires to define ALLOW_SITRACER # define ALLOW_SITRACER -# define ALLOW_SITRACER_ADVCAP +#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-- Allow sea-ice dynamic code. These options are 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). 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. +C discretization should still work but it is not longer actively +C tested and supported. Define this flag to compile it. It cannot be +C defined together with SEAICE_CGRID +#undef SEAICE_BGRID_DYNAMICS + +C-- The following flag should always be set in order to use C the +C-- operational C-grid discretization. #define SEAICE_CGRID -C-- Only for the C-grid version it is possible to #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 +# define 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- 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 + +CMLC or the modified coulombic rheology +CMLC allow the use of the Mohr Coulomb rheology (runtime flag +CMLC SEAICEuseFULLMC) as defined in (Ip 1991) /!\ This is known +CMLC to give unstable results, use with caution +CML# define SEAICE_ALLOW_FULLMC + +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 @@ -122,56 +171,108 @@ 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. # define 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-- By default for B-grid dynamics solver wind stress under sea-ice is + +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS +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- 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- 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-- When set limit the Ice-Loading to mass of 1/5 of Surface ocean grid-box +#endif /* SEAICE_BGRID_DYNAMICS */ + +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- 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- 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 +C This options excludes more complex physics such +C as sublimation, ITD, and frazil. #define 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- Do not compile/use seaice-related obcs code when using obcs. +#undef DISABLE_SEAICE_OBCS + C-- Enable free drift code #define SEAICE_ALLOW_FREEDRIFT C-- pkg/seaice cost functions compile flags -c >>> Sea-ice volume (requires pkg/cost) +C- Sea-ice volume (requires pkg/cost) #define 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 */ - -CEH3 ;;; Local Variables: *** -CEH3 ;;; mode:fortran *** -CEH3 ;;; End: *** diff --git a/verification/seaice_itd/code/SEAICE_OPTIONS.h b/verification/seaice_itd/code/SEAICE_OPTIONS.h index f9697ddec7..dd6c75f73e 100644 --- a/verification/seaice_itd/code/SEAICE_OPTIONS.h +++ b/verification/seaice_itd/code/SEAICE_OPTIONS.h @@ -1,3 +1,8 @@ +#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. @@ -6,24 +11,12 @@ C | Use this file for selecting options within the sea ice C | package. C *==========================================================* -#ifndef SEAICE_OPTIONS_H -#define SEAICE_OPTIONS_H -#include "PACKAGES_CONFIG.h" -#include "CPP_OPTIONS.h" - #ifdef ALLOW_SEAICE -C Package-specific Options & Macros go here +C--- Package-specific Options & Macros go here C-- Write "text-plots" of certain fields in STDOUT for debugging. #undef SEAICE_DEBUG -C-- Allow sea-ice dynamic code. -C This option is provided to allow use of TAMC -C on the thermodynamics component of the code only. -C Sea-ice dynamics can also be turned off at runtime -C using variable SEAICEuseDYNAMICS. -#define SEAICE_ALLOW_DYNAMICS - 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 @@ -55,68 +48,116 @@ 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-- Tracers of ice and/or ice cover. -#undef ALLOW_SITRACER -#ifdef ALLOW_SITRACER -C-- To try avoid 'spontaneous generation' of tracer maxima by advdiff. -# define ALLOW_SITRACER_ADVCAP -#endif - -C-- Enable grease ice parameterization -C The grease ice parameterization delays formation of solid -C sea ice from frazil ice by a time constant and provides a -C dynamic calculation of the initial solid sea ice thickness -C HO as a function of winds, currents and available grease ice -C volume. Grease ice does not significantly reduce heat loss -C from the ocean in winter and area covered by grease is thus -C handled like open water. -C (For details see Smedsrud and Martin, 2014, Ann.Glac.) +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 -C to yield grease ice volume. Additionally, the actual grease ice -C layer thickness (diagnostic SIgrsLT) can be saved. +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-- grease ice uses SItracer: + +C-- Tracers of ice and/or ice cover. #ifdef SEAICE_GREASE +C SEAICE_GREASE code requires to define ALLOW_SITRACER # define ALLOW_SITRACER -# define ALLOW_SITRACER_ADVCAP +#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-- Allow sea-ice dynamic code. These options are 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). 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. +C discretization should still work but it is not longer actively +C tested and supported. Define this flag to compile it. It cannot be +C defined together with SEAICE_CGRID +#undef SEAICE_BGRID_DYNAMICS + +C-- The following flag should always be set in order to use C the +C-- operational C-grid discretization. #define SEAICE_CGRID -C-- Only for the C-grid version it is possible to #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 # define SEAICE_ALLOW_JFNK + C enable Krylov code by defining the following flag # define SEAICE_ALLOW_KRYLOV -C enable this flag to reproduce old verification results for JFNK + +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 # define SEAICE_PRECOND_EXTRA_EXCHANGE + C enable LSR to use global (multi-tile) tri-diagonal solver # define 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- 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 # define 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 @@ -124,54 +165,108 @@ 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-- By default for B-grid dynamics solver wind stress under sea-ice is + +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS +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- 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- 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-- When set limit the Ice-Loading to mass of 1/5 of Surface ocean grid-box +#endif /* SEAICE_BGRID_DYNAMICS */ + +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- 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- 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 +C This options excludes more complex physics such +C as sublimation, ITD, and frazil. +#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- Do not compile/use seaice-related obcs code when using obcs. +#undef DISABLE_SEAICE_OBCS + C-- Enable free drift code #define SEAICE_ALLOW_FREEDRIFT C-- pkg/seaice cost functions compile flags -c >>> Sea-ice volume (requires pkg/cost) +C- Sea-ice volume (requires pkg/cost) #undef ALLOW_COST_ICE -c >>> Sea-ice misfit to obs (requires pkg/cost and ecco) -#undef ALLOW_SEAICE_COST_SMR_AREA +#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 */ - -CEH3 ;;; Local Variables: *** -CEH3 ;;; mode:fortran *** -CEH3 ;;; End: *** diff --git a/verification/seaice_obcs/code/SEAICE_OPTIONS.h b/verification/seaice_obcs/code/SEAICE_OPTIONS.h index 2ca9e5012e..83350b5875 100644 --- a/verification/seaice_obcs/code/SEAICE_OPTIONS.h +++ b/verification/seaice_obcs/code/SEAICE_OPTIONS.h @@ -1,3 +1,8 @@ +#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. @@ -6,24 +11,12 @@ C | Use this file for selecting options within the sea ice C | package. C *==========================================================* -#ifndef SEAICE_OPTIONS_H -#define SEAICE_OPTIONS_H -#include "PACKAGES_CONFIG.h" -#include "CPP_OPTIONS.h" - #ifdef ALLOW_SEAICE -C Package-specific Options & Macros go here +C--- Package-specific Options & Macros go here C-- Write "text-plots" of certain fields in STDOUT for debugging. #undef SEAICE_DEBUG -C-- Allow sea-ice dynamic code. -C This option is provided to allow use of TAMC -C on the thermodynamics component of the code only. -C Sea-ice dynamics can also be turned off at runtime -C using variable SEAICEuseDYNAMICS. -#define SEAICE_ALLOW_DYNAMICS - 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 @@ -55,66 +48,116 @@ C seaice with a fraction (=SEAICE_saltFrac) of freezing seawater salinity. C- Note: SItracer also offers an alternative way to handle variable salinity. #define SEAICE_VARIABLE_SALINITY -C-- Tracers of ice and/or ice cover. -#undef ALLOW_SITRACER -#ifdef ALLOW_SITRACER -C-- To try avoid 'spontaneous generation' of tracer maxima by advdiff. -# define ALLOW_SITRACER_ADVCAP -#endif - -C-- Enable grease ice parameterization -C The grease ice parameterization delays formation of solid -C sea ice from frazil ice by a time constant and provides a -C dynamic calculation of the initial solid sea ice thickness -C HO as a function of winds, currents and available grease ice -C volume. Grease ice does not significantly reduce heat loss -C from the ocean in winter and area covered by grease is thus -C handled like open water. -C (For details see Smedsrud and Martin, 2014, Ann.Glac.) +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 -C to yield grease ice volume. Additionally, the actual grease ice -C layer thickness (diagnostic SIgrsLT) can be saved. +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-- grease ice uses SItracer: + +C-- Tracers of ice and/or ice cover. #ifdef SEAICE_GREASE +C SEAICE_GREASE code requires to define ALLOW_SITRACER +# define ALLOW_SITRACER +#else # define ALLOW_SITRACER -# define ALLOW_SITRACER_ADVCAP #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-- Allow sea-ice dynamic code. These options are 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). 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. +C discretization should still work but it is not longer actively +C tested and supported. Define this flag to compile it. It cannot be +C defined together with SEAICE_CGRID +#undef SEAICE_BGRID_DYNAMICS + +C-- The following flag should always be set in order to use C the +C-- operational C-grid discretization. #define SEAICE_CGRID -C-- Only for the C-grid version it is possible to #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 # define SEAICE_ALLOW_JFNK + C enable Krylov code by defining the following flag # define 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- 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 @@ -122,54 +165,108 @@ 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-- By default for B-grid dynamics solver wind stress under sea-ice is + +#endif /* SEAICE_CGRID */ + +#ifdef SEAICE_BGRID_DYNAMICS +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- 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- 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-- When set limit the Ice-Loading to mass of 1/5 of Surface ocean grid-box +#endif /* SEAICE_BGRID_DYNAMICS */ + +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- 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- 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 +C This options excludes more complex physics such +C as sublimation, ITD, and frazil. +#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- Do not compile/use seaice-related obcs code when using obcs. +#undef DISABLE_SEAICE_OBCS + C-- Enable free drift code #define SEAICE_ALLOW_FREEDRIFT C-- pkg/seaice cost functions compile flags -c >>> Sea-ice volume (requires pkg/cost) +C- Sea-ice volume (requires pkg/cost) #undef ALLOW_COST_ICE -c >>> Sea-ice misfit to obs (requires pkg/cost and ecco) -#undef ALLOW_SEAICE_COST_SMR_AREA +#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 */ - -CEH3 ;;; Local Variables: *** -CEH3 ;;; mode:fortran *** -CEH3 ;;; End: ***