bkgf / romsiceshelf Goto Github PK

View Code? Open in Web Editor NEW

8.0 8.0 2.0 60.64 MB

romsIceShelf

Makefile 0.67% C++ 6.01% Fortran 77.27% C 15.30% TeX 0.05% Shell 0.48% Awk 0.01% Perl 0.21% SourcePawn 0.01%

romsiceshelf's People

Contributors

Stargazers

Watchers

Forkers

dgwyther mmainak

romsiceshelf's Issues

Hz does not have a type

Following the previous issue (ANA_SEAICE), Hz does not seem to have a type.
Modifications need to be considered in Nonlinear/Iceshelf/iceshelf_vbc.F.

I have committed a fix for this on my brach called 'mertz', and commit called 'Hz definition'. (See the README.txt)

Periodic boundary conditions obsolete CPP tags

The periodic boundary condition CPP options and logicals as defined in ROMS/Nonlinear/IceShelf/iceshelf_vbc.F at ~line 178 are obsolete. Periodic boundary conditions are now defined as on/off in the *.in file. As a result, these lines throw up errors.

To remove this error, the lines:

 ifdef DISTRIBUTE
#  ifdef EW_PERIODIC
      logical :: EWperiodic=.TRUE.
#  else
      logical :: EWperiodic=.FALSE.
#  endif
#  ifdef NS_PERIODIC
      logical :: NSperiodic=.TRUE.
#  else
      logical :: NSperiodic=.FALSE.
#  endif
# endif

should be removed/commented out.

zice "is not a field name that is defined in the encompassing structure"

Error while compiling romsIceShelf:

step2d.f90(77): error #6460: This is not a field name that is defined in the encompassing structure.   [ZICE]
     &                  GRID(ng) % zice,                                &
-----------------------------------^
step2d.f90(112): error #6404: This name does not have a type, and must have an explicit type.   [ZICE]
     &                        zice,                                     &
------------------------------^

Upon investigation, this line originates from ROMS/Nonlinear/step2d_LF_AM3.h, line 63:

# ifdef SOLVE3D && defined ICESHELF
     &                  GRID(ng) % zice,                                &

Which is inserted into step2d.F when the CPP does its magic. I was running the upwelling test case, /without/ #define ICESHELF in upwelling.h. Note that when I add this option to upwelling.h, the issue is resolved. (However, other issues involving iceshelf_vbc.F are raised). So, before the iceshelf mods can be merged back into the ROMS trunk, step2d_LF_AM3.h will need to be modified to allow a non-iceshelf model run.

Bug - Basal T and S; bad code

In the calculation of the iceshelf vertical boundary conditions, in the file ROMS/Nonlinear/IceShelf/iceshelf_vbc.F, the calculation of the basal temperature contains an error: mflag*(cp_i/cp_w)*m(i,j)+). The '+)' is either an unfinished equation or a mistyped line.

# ifdef ICESHELF_VBC_HFPROPTOM
! Calculates density RHOW in the boundary layer, interface
! pressure PG[dbar], and solving a quadratic equation for the 
! interface salinity (SB) to determine the melTemp_insitug/freezingrate
! (seta).
...
# else
...
! Calculate basal temperature and salinity:
           Tb = (gammaT*Tm+mflag*(cp_i/cp_w)*m(i,j)*Ti-(L/cp_w)*m(i,j)) &
              / (gammaT + mflag*(cp_i/cp_w)*m(i,j)+)
           Sb = (Tb - b - c*zice(i,j))/a
...

('...' indicate omitted lines of code)

ROMS/Nonlinear/zetabc.F

In this file there is a bug in Chapman boundary conditions for zeta (both implicit and explicit, and along all four boundaries). When ice shelves are turned on, the zeta term is missing in the calculation of cff1. The fixed code is in the repository https://github.com/bkgf/ROMSIce in the folder ROMSIceShelf.

BSTRESS, ICESTRESS?

BSTRESS and ICESTRESS options seem to increase stability. Should be pushed into master branch.

Bug - Ustar calculated from sustr and svstr, before sustr and svstr calculated.

I noticed that the friction velocity, Ustar, is calculated from the ice shelf surface stresses, sustr and svstr, BEFORE sustr and svstr are calculated. This is in ROMS/Nonlinear/IceShelf/iceshelf_vbc.F. This seems to suggest that the first calculation of Ustar could be calculated with the sustr and svstr which are calculated analytically from surface wind stress across the open ocean.

Calculation of adjusted Prandtl and Schmidt numbers

In iceshelf_vbc.F, the adjustment to the Prandtl and Schmidt numbers,

Pradj = 12.5_r8*(Pr**(2/3))
Scadj = 12.5_r8*(Sc**(2/3))

Uses the wrong number declaration. As it is, they are declared as integers (by default), and so the adjusted Prandtl and Schmidt numbers evaluate to (both) 12.5. (As 2/3 evaluates and truncates to 0).

Consequently, this line should read:

Pradj = 12.5_r8*(Pr**(2.0_r8/3.0_r8))
Scadj = 12.5_r8*(Sc**(2.0_r8/3.0_r8))

(Or by using some other way of declaring them as real. i.e. 2d0)

bad ice density gradient in prsgrd32.h

In ROMS/Nonlinear/prsgrd32.h:
Change

#ifdef ICESHELF
      real(r8), parameter :: drhodz = 0.0_r8 !0.00478_r8
#endif

#ifdef ICESHELF
      real(r8), parameter :: drhodz = 0.00478_r8
#endif

Issue - Very slow heat exchange

As it is, heat exchange will be very slow in the idealised case. This is due to several issues. The first two are raised in this issue report.

[ ] In ROMS/Nonlinear/Iceshelf/iceshelf_vbc.F, if ustar is small, it is set to a value small;

           if(ustar.lt.small) ustar = small

However, the value of small is by default (1 x 10^-10) so small, that heat exchange due to ustar will be SO low, that buoyancy driven currents will be weak and ustar will never increase above the value of small. It is better to put the lower bound of small higher - say 1 x 10^-4. This can be affected in iceshelf_mod.h. (see my commit below)

[ ] Another issue is that in iceshelf_vbc.F, logic dictates that if ustar is smaller than small, it is set to small. As stated above, if this is the case, heat exchange will be very slow, and ustar will take a LONG time to increase above small. Consequently, turb will always = 0, and heat exchange will be further discouraged. This suppression of heat exchange can be removed by changing the logic gate to capture the case where ustar = small:

            IF (ustar.gt.small.and.ABS(f(i,j)).gt.1.0E-8) THEN
              turb = 2.5*LOG(5300.0_r8*ustar*ustar/ABS(f(i,j))) + 7.12
            ELSE
              turb = 0.0_r8
            END IF

So here, if ustar < small, ustar = ustar, then, turbulence will be set to 0, reinforcing the suppression of heat exchange. If we instead change the greater than to be greater than or equal to, turbulence will be included.

My edits which fix this are shown in: dgwyther@290264e

write out melt rate after restart

If you restart from a previous solution you can get this error:
WRT_HIS - error while writing variable: m
After deactivating all ice shelf write out options, the restart works fine.

Example selection of time step and write out times to reproduce the error:

Input/Output parameters.

       NRREC == -1
   LcycleRST == T
        NRST == 10
...

! Output history, average, diagnostic files parameters.

     LDEFOUT == T
        NHIS == 9 
     NDEFHIS == 315360

Bug in use of surface fluxes.

ROMSIceShelf/ROMS/Nonlinear/IceShelf/iceshelf_vbc.F

Lines 561 to 564 in b95499a

    
           # else 
        
                     DO itrc=1,NT(ng) 
        
                        stflx(i,j,itrc)=0.0_r8 
        
                     END DO

These lines, when not using the flag ANA_SEAICE, will result in ALL surface heat and water forcing to be set to 0. These four lines should be commented out.

Bug - Calculation of quadratic drag stresses; bad code.

There is an error in the calculation of the drag stresses, sustr and svstr, in the quadratic parameterisation. This is a coding error, where cff2 is used instead of cff1. The code has been corrected in my branch.

Missing MASKING option

Probably need to have this in ROMS/Include/amery.h:

#define MASKING

ROMSIceShelf/ROMS/External/ocean_icetest.in

Instead of
MyAppCPP = WEDDELL

it should be
MyAppCPP = ICETEST

update .in files for new format

Has anyone thought about updating the iceshelf specific .in files to the current format? I can see that ocean_icetest.in has been updated, but ocean_amery.in hasn't, nor has ocean_waomX.in.
@mdinniman I can update and make a push request.

Possible issue - CPP label not recognised

I'm trying to run the ICETEST scenario. I have added the CPP definition #define ICESHELF_VBC_HFPROPTOM to icetest.h (it is included below).

I am doing this as I want C to remove the block of code that includes the bug I previously reported. As this bug is present in the #else block, and I defined ICESHELF_VBC_HFPROPTOM, then the buggy code block should not be compiled into the file iceshelf_vbc.f90.

However, even though I've defined this option, the CPP is still included the #else block, indicating that it's not recognising the CPP label ICESHELF_VBC_HFPROPTOM.

From iceshelf_vbc.F

# ifdef ICESHELF_VBC_HFPROPTOM

! Calculates density RHOW in the boundary layer, interface
! pressure PG[dbar], and solving a quadratic equation for the 
! interface salinity (SB) to determine the melTemp_insitug/freezingrate
! (seta).
          rho_r= rho_i/(1000.0_r8+rho(i,j,N(ng)))
          Ep1  = cp_w*gammaT
          Ep2  = cp_i*gammaS
          Ep3  = L*gammaS
          Ep31 = -rho_r*cp_i*dt_i/zice(i,j)
          Ep4  = b+c*z_w(i,j,N(ng))
          Ep5  = gammaS/rho_r
!
! Temperature gradient within the ice ONLY changes with melTemp_insitug
          TFb = a*Sm+Ep4
          IF (Tm.lt.TFb) THEN
            Ex1 = a*(Ep1+Ep31)
            Ex2 = Ep1*(Tm-Ep4)+Ep3+Ep31*(Ti-Ep4)
            Ex3 = Ep3*Sm
            Ex6 = 0.5
          ELSE
! negative heat flux term in the ice (due to -kappa/D)
            Ex1 = a*(Ep1-Ep2)
            Ex2 = Ep1*(Ep4-Tm)+Ep2*(Ti+a*Sm-Ep4)-Ep3
            Ex3 = Sm*(Ep2*(Ep4-Ti)+Ep3)
            Ex6 = -0.5
          endif
          Ex4 = Ex2/Ex1
          Ex5 = Ex3/Ex1
          Sr1 = 0.25*Ex4*Ex4-Ex5
          Sr2 = Ex6*Ex4
          Sf1 = Sr2+sqrt(Sr1)
          Tf1 = a*Sf1+Ep4
          Sf2 = Sr2-sqrt(Sr1)
          Tf2 = a*Sf2+Ep4
!
! Salinities < 0 psu are NOT defined, i.e....
         IF (Sf1.gt.0.) THEN
           Tb = Tf1
           Sb = Sf1
         ELSE
           Tb = Tf2
           Sb = Sf2
         endif
!
         m(i,j) = Ep5*(1.0-Sm/Sb)
!         m(i,j)= Seta * dtfast
!
! Calculates fluxes:
         stflx(i,j,itemp)=(gammaT-rho_r*m(i,j))*(Tb-Tm)
         stflx(i,j,isalt)=(gammaS-rho_r*m(i,j))*(Sb-Sm)

# else
            TFb = a*Sm + b + c*zice(i,j)
            mflag= 0.5*(1 + SIGN(1.0_r8,Tm-TFb))
            TFi = (1 - mflag)*a*Si + b + c*zice(i,j)

! Calculate coefficents in quadratic to be solved:
           cff1 = L/cp_w + mflag*(cp_i/cp_w)*(TFi-Ti)
           cff2 = gammaS*(L/cp_w + mflag*(cp_i/cp_w)*(TFb-Ti))          &
     &          + gammaT*(TFi-Tm)
           cff3 = gammaS*gammaT*(TFb - Tm)

! Calculate melt rate:
           m(i,j) = -(cff2 - SQRT(cff2*cff2 - 4*cff1*cff3))/(2*cff1)

! Calculate basal temperature and salinity:
           Tb = (gammaT*Tm+mflag*(cp_i/cp_w)*m(i,j)*Ti-(L/cp_w)*m(i,j)) &
              / (gammaT + mflag*(cp_i/cp_w)*m(i,j)+)
           Sb = (Tb - b - c*zice(i,j))/a

           stflx(i,j,itemp)=(gammaT+rhoi_on_rho0*m(i,j))*(Tb-Tm)
           stflx(i,j,isalt)=(gammaS+rhoi_on_rho0*m(i,j))*(Sb-Sm)
!           write(6,*) gammaT,gammaS,Tb,Sb,Tm,Sm,rhoi_on_rho0
# endif

my icetest.h:

** Copyright (c) 2002-2008 The ROMS/TOMS Group                               **
**   Licensed under a MIT/X style license                                    **
**   See License_ROMS.txt                                                    **
*******************************************************************************
**
** Options for Simplified Ice Shelf Ocean Cavity Model Test.
**
** Application flag:   ICETEST
** Input script:       ocean_icetest.in
*/
#define UV_ADV
#define DJ_GRADPS
#define UV_COR
#define UV_VIS2
#define UV_QDRAG
#define MIX_GEO_UV
#undef MIX_S_UV
#define TS_C4HADVECTION
#define TS_C4VADVECTION
#define TS_DIF2
#define MIX_GEO_TS
#undef  MIX_S_TS
#define SOLVE3D
#define SALINITY
#define NONLIN_EOS
#define CURVGRID
#define SPHERICAL
#define SPLINES
#define ICESHELF
#undef ICESHELF_MORPH
#define ICESHELF_VBC_HFPROPTOM
#undef  AVERAGES
#undef ATM_PRESS
#undef ANA_PAIR
#define PERFECT_RESTART
#define ANA_GRID
#define ANA_INITIAL
#define ANA_SMFLUX
#define ANA_STFLUX
#define ANA_SSFLUX
#define ANA_BSFLUX
#define ANA_BTFLUX
#define ANA_SRFLUX

/* Define SET_VBC.F for open ocean boundary layer. Can be one of:
* * ANA_SEAICE
*  Note that both undef will set surface fluf of salt and temp to zero*/
#define ANA_SEAICE
/* Define SET_VBC.F for ice-ocean Thermodynamics. Can be one of:
*  * VBC_ICE_2EQN
*  * VBC_ICE_3EQN       
*  * Note that both undef will set surface fluf of salt and temp to zero */
#undef ICESHELF_2EQN_VBC
#define ICESHELF_3EQN_VBC
#undef ICESHELF_TEOS10


#undef  ANA_VMIX
#undef MY25_MIXING
#define LMD_MIXING

#ifdef MY25_MIXING
#define N2S2_HORAVG
#define K_C4ADVECTION
#endif

#ifdef  LMD_MIXING
#define LMD_CONVEC 
#undef LMD_DDMIX 
#undef LMD_RIMIX
#undef LMD_SKPP
#endif

#undef SEDIMENT
#ifdef SEDIMENT
#define ANA_SEDIMENT
#define SED_MORPH
/*#define SED_DENS
#define SUSPLOAD
#define ANA_BPFLUX
#define ANA_SPFLUX*/
#endif

/* Test groundwater fluxes */
#undef UV_PSOURCE
#undef ANA_PSOURCE
#undef TS_PSOURCE

When I run, using makefile with ICETEST as the application, I get the following at the command line:

cd Build; /apps/sles11sp2/intel/composer_xe_2013.1.117/bin/intel64/ifort -c  -g -check bounds -traceback -check uninit -warn interfaces,nouncalled -gen-interfaces iceshelf_vbc.f90
iceshelf_vbc.f90(269): error #5082: Syntax error, found ')' when expecting one of: ( <IDENTIFIER> <CHAR_CON_KIND_PARAM> <CHAR_NAM_KIND_PARAM> <CHARACTER_CONSTANT> <INTEGER_CONSTANT> ...
              / (gammaT + mflag*(cp_i/cp_w)*m(i,j)+)
---------------------------------------------------^
compilation aborted for iceshelf_vbc.f90 (code 1)
make: *** [Build/iceshelf_vbc.o] Error 1

And here is the relevant lines of iceshelf_vbc.f90:

          IF (zice(i,j).ne.0.0_r8) THEN
! Get the insitu salinity and temperature
         CALL potit(Sm,t(i,j,N(ng),nrhs,itemp),-zice(i,j),0.0_r8,Tm,i,j)
!    Calculate exchange coefficients. Note that this assumes log profile 
            ustar = SQRT(SQRT((0.5_r8*(sustr(i,j)+sustr(i+1,j)))**2+    &
     &                        (0.5_r8*(svstr(i,j)+svstr(i,j+1)))**2))
           if(ustar.lt.small) ustar = small
! Uses simplified version of McPhee 1987
            IF (ustar.gt.small.and.ABS(f(i,j)).gt.1.0E-8) THEN
              turb = 2.5*LOG(5300.0_r8*ustar*ustar/ABS(f(i,j))) + 7.12
            ELSE
              turb = 0.0_r8
            END IF
            gammaT = ustar/(turb + Pradj - 6.0_r8)
            gammaS = ustar/(turb + Scadj - 6.0_r8)
            TFb = a*Sm + b + c*zice(i,j)
            mflag= 0.5*(1 + SIGN(1.0_r8,Tm-TFb))
            TFi = (1 - mflag)*a*Si + b + c*zice(i,j)
! Calculate coefficents in quadratic to be solved:
           cff1 = L/cp_w + mflag*(cp_i/cp_w)*(TFi-Ti)
           cff2 = gammaS*(L/cp_w + mflag*(cp_i/cp_w)*(TFb-Ti))          &
     &          + gammaT*(TFi-Tm)
           cff3 = gammaS*gammaT*(TFb - Tm)
! Calculate melt rate:
           m(i,j) = -(cff2 - SQRT(cff2*cff2 - 4*cff1*cff3))/(2*cff1)
! Calculate basal temperature and salinity:
           Tb = (gammaT*Tm+mflag*(cp_i/cp_w)*m(i,j)*Ti-(L/cp_w)*m(i,j)) &
              / (gammaT + mflag*(cp_i/cp_w)*m(i,j)+)
           Sb = (Tb - b - c*zice(i,j))/a
           stflx(i,j,itemp)=(gammaT+rhoi_on_rho0*m(i,j))*(Tb-Tm)
           stflx(i,j,isalt)=(gammaS+rhoi_on_rho0*m(i,j))*(Sb-Sm)
!           write(6,*) gammaT,gammaS,Tb,Sb,Tm,Sm,rhoi_on_rho0
          END IF

stflx is not calculated with ANA_SEAICE for the open ocean

As it is, there is no surface fluxes calculated over the open ocean. In this issue report I list the modifications I have made on my edits.

First modification (in Nonlinear/Iceshelf/iceshelf_mod.F).
As it is now the variables 'gamma', 'refSalt, and 'temp_f' used for ICESHELF_3EQN_VBC are defined only if ANA_SEAICE is defined.

#   ifdef ICESHELF_3EQN_VBC
      real(r8), parameter :: gamma = 0.0001_r8
      real(r8), parameter :: refSalt = 34.4_r8
      real(r8) :: temp_f
#   endif

This needs to be moved earlier in the script where the other variables for ICESHELF_3EQN_VBC are defined.

Second modifications (Nonlinear/IceShelf/iceshelf_vbc.F).
As it is now, stflx over the open ocean is not calculated. Line 266 a IF statement for zice =~ 0 (where there is an ice shelf) starts and is followed by another IF statement (l.267) which lists all the different ways of calculating the fluxes at the ice/ocean interphase. However, this statement finishes line 423, after the second second block of the IF statement which consider the case of the open ocean (line 406 -- ELSE).

         IF (zice(i,j).ne.0.0_r8) THEN
#if defined ICESHELF_2EQN_VBC
          TFb = a*Sm+b+c*zice(i,j)
…

          ELSE
            m(i,j) = 0.0_r8
# if defined ANA_SEAICE
            stflx(i,j,itemp)=Hz(i,j,N(ng))*                             &
                             (sfcTemp-t(i,j,N(ng),nrhs,itemp))/trelax
            stflx(i,j,isalt)=Hz(i,j,N(ng))*                             &
                             (sfcSalt-t(i,j,N(ng),nrhs,isalt))/trelax
# else
          DO itrc=1,NT(ng)
             stflx(i,j,itrc)=0.0_r8
          END DO
# endif

!          END IF
!        END DO
!      END DO

#endif

          END IF

This means stflx from the ANA_SEAICE is not calculated for the open ocean. #endif currently at line 423 needs to be moved before the ELSE at line 406.

I have committed a fix for this. See my branch called 'mertz' and then 'ANA_SEAICE' commit. Details are given in the README.txt file.