shelfice_remeshing/code/shelfice_massmin.F

C $Header: /u/gcmpack/MITgcm_contrib/verification_other/shelfice_remeshing/code/shelfice_massmin.F,v 1.4 2017/02/10 18:52:46 dgoldberg Exp $
C $Name:  $

C KS_dens------

#include "SHELFICE_OPTIONS.h"

CBOP
      SUBROUTINE SHELFICE_MASSMIN( R_min, massMin, myThid )
C     *============================================================*
C     | SUBROUTINE SHELFICE_MASSMIN
C     | o Routine to determine the minimum mass based on MCWT
C     *============================================================*
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SHELFICE.h"
#include "DYNVARS.h"
#ifdef ALLOW_COST
# include "SHELFICE_COST.h"
#endif /* ALLOW_COST */

C     === Routine arguments ===
C     myThid -  Number of this instance 
C     
C     This subroutine calculates the minimum mass to allow MWCT space
C     under the grounded ice.

      _RL massMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL R_min(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      INTEGER bi
      INTEGER bj
      INTEGER myThid

CEndOfInterface

#ifdef ALLOW_SHELFICE
C     === Local variables ===
C     i,j,bi,bj - Loop counters
      INTEGER i, j, k, k2
      _RL drho (Nr)
      _RL phiHydF(Nr+1)
      _RL phiHydC(Nr)
      _RL vert_rho(Nr)
      _RL drLoc, pLoad_min, sum_hFac
      _RL SEALEVEL, ETA
CEOP

C Initialize variables
      R_min = 0 _d 0
      drLoc = 0 _d 0
      pLoad_min = 0 _d 0
      DO k=1,Nr
       vert_rho(k) = 0 _d 0
      ENDDO
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          massMin(i,j,bi,bj) = 0 _d 0
         ENDDO
        ENDDO
       ENDDO
      ENDDO

C First need vertical profile of density anomaly;
      CALL SHELFICE_VERT_DENS( vert_rho, myThid )
C Get the average open ocean ssh
      SEALEVEL = 0. _d 0
      ETA = 0. _d 0
      CALL SHELFICE_SEA_LEVEL_AVG( SEALEVEL,ETA, myThid )
C calcaulte hydrostatic pressure based on density profile
      DO k = 1,Nr
       drho(k) = vert_rho(k)
        IF (k.EQ.1) THEN
          phiHydF(k) = 0.
C          phiHydF(k) = SEALEVEL*gravity*drho(k)*recip_rhoConst
        ENDIF
          phiHydC(k)=phiHydF(k) + halfRL*drF(k)*gravity*drho(k)
     &                                        *recip_rhoConst
          phiHydF(k+1)=phiHydC(k) + halfRL*drF(k)*gravity*drho(k)
     &                                        *recip_rhoConst
      ENDDO

C calculate the level of minimum pressure and from that massMin
C based on R_MCWT
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1,sNy
          DO i=1,sNx
           R_min(i,j,bi,bj)=R_Low(i,j,bi,bj)+R_MWCT(i,j,bi,bj)
           DO k =1,Nr
            IF ((R_min(i,j,bi,bj) .GT. rF(k+1)) .AND.
     &          (R_min(i,j,bi,bj) .LE. rF(k)))  THEN
              drLoc = (rF(k)-R_min(i,j,bi,bj))/drF(k)

             pLoad_min = (phiHydF(k) + drLoc*
     &           (phiHydF(k+1)-phiHydF(k)))*rhoConst
     &           + gravity*rhoConst*(drLoc*drF(k))

             massMin(i,j,bi,bj) = pLoad_min*recip_gravity - 
     &             rF(k)*rhoConst
            ENDIF
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDDO


#endif /* ALLOW_SHELFICE */

      RETURN
      END
1	dgoldberg	1.5	C $Header: /u/gcmpack/MITgcm_contrib/verification_other/shelfice_remeshing/code/shelfice_massmin.F,v 1.4 2017/02/10 18:52:46 dgoldberg Exp $
2	ksnow	1.1	C $Name: $
3
4			C KS_dens------
5
6			#include "SHELFICE_OPTIONS.h"
7
8			CBOP
9	dgoldberg	1.4	SUBROUTINE SHELFICE_MASSMIN( R_min, massMin, myThid )
10	ksnow	1.1	C ============================================================
11			C \| SUBROUTINE SHELFICE_MASSMIN
12			C \| o Routine to determine the minimum mass based on MCWT
13			C ============================================================
14			IMPLICIT NONE
15
16			C === Global variables ===
17			#include "SIZE.h"
18			#include "EEPARAMS.h"
19			#include "PARAMS.h"
20			#include "GRID.h"
21			#include "SHELFICE.h"
22			#include "DYNVARS.h"
23			#ifdef ALLOW_COST
24			# include "SHELFICE_COST.h"
25			#endif /* ALLOW_COST */
26
27			C === Routine arguments ===
28			C myThid - Number of this instance
29			C
30			C This subroutine calculates the minimum mass to allow MWCT space
31			C under the grounded ice.
32
33			_RL massMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
34			_RL R_min(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
35	dgoldberg	1.3	INTEGER bi
36			INTEGER bj
37	ksnow	1.1	INTEGER myThid
38
39			CEndOfInterface
40
41			#ifdef ALLOW_SHELFICE
42			C === Local variables ===
43			C i,j,bi,bj - Loop counters
44	dgoldberg	1.3	INTEGER i, j, k, k2
45	ksnow	1.1	_RL drho (Nr)
46			_RL phiHydF(Nr+1)
47			_RL phiHydC(Nr)
48			_RL vert_rho(Nr)
49			_RL drLoc, pLoad_min, sum_hFac
50	ksnow	1.2	_RL SEALEVEL, ETA
51	ksnow	1.1	CEOP
52
53			C Initialize variables
54			R_min = 0 _d 0
55			drLoc = 0 _d 0
56			pLoad_min = 0 _d 0
57			DO k=1,Nr
58			vert_rho(k) = 0 _d 0
59			ENDDO
60	dgoldberg	1.5	DO bj = myByLo(myThid), myByHi(myThid)
61			DO bi = myBxLo(myThid), myBxHi(myThid)
62	ksnow	1.1	DO j=1-OLy,sNy+OLy
63			DO i=1-OLx,sNx+OLx
64			massMin(i,j,bi,bj) = 0 _d 0
65			ENDDO
66			ENDDO
67	dgoldberg	1.5	ENDDO
68			ENDDO
69	ksnow	1.1
70			C First need vertical profile of density anomaly;
71			CALL SHELFICE_VERT_DENS( vert_rho, myThid )
72	ksnow	1.2	C Get the average open ocean ssh
73			SEALEVEL = 0. _d 0
74			ETA = 0. _d 0
75			CALL SHELFICE_SEA_LEVEL_AVG( SEALEVEL,ETA, myThid )
76	ksnow	1.1	C calcaulte hydrostatic pressure based on density profile
77			DO k = 1,Nr
78			drho(k) = vert_rho(k)
79			IF (k.EQ.1) THEN
80			phiHydF(k) = 0.
81	ksnow	1.2	C phiHydF(k) = SEALEVELgravitydrho(k)*recip_rhoConst
82	ksnow	1.1	ENDIF
83			phiHydC(k)=phiHydF(k) + halfRLdrF(k)gravity*drho(k)
84			& *recip_rhoConst
85			phiHydF(k+1)=phiHydC(k) + halfRLdrF(k)gravity*drho(k)
86			& *recip_rhoConst
87	dgoldberg	1.3	ENDDO
88	ksnow	1.1
89			C calculate the level of minimum pressure and from that massMin
90			C based on R_MCWT
91	dgoldberg	1.4	DO bj = myByLo(myThid), myByHi(myThid)
92			DO bi = myBxLo(myThid), myBxHi(myThid)
93	dgoldberg	1.3	DO j=1,sNy
94			DO i=1,sNx
95			R_min(i,j,bi,bj)=R_Low(i,j,bi,bj)+R_MWCT(i,j,bi,bj)
96	ksnow	1.1	DO k =1,Nr
97			IF ((R_min(i,j,bi,bj) .GT. rF(k+1)) .AND.
98			& (R_min(i,j,bi,bj) .LE. rF(k))) THEN
99	ksnow	1.2	drLoc = (rF(k)-R_min(i,j,bi,bj))/drF(k)
100	ksnow	1.1
101			pLoad_min = (phiHydF(k) + drLoc*
102			& (phiHydF(k+1)-phiHydF(k)))*rhoConst
103	ksnow	1.2	& + gravityrhoConst(drLoc*drF(k))
104
105	ksnow	1.1	massMin(i,j,bi,bj) = pLoad_min*recip_gravity -
106	ksnow	1.2	& rF(k)*rhoConst
107	ksnow	1.1	ENDIF
108			ENDDO
109			ENDDO
110			ENDDO
111	dgoldberg	1.4	ENDDO
112			ENDDO
113	dgoldberg	1.3
114	ksnow	1.1
115			#endif /* ALLOW_SHELFICE */
116
117			RETURN
118			END