Merge pull request #10 from dfer22/bugfix_thebug

Code for MSE coordinate (first try with Github!)

Author: jm-c

doc/tag-index

@@ -1,6 +1,7 @@
     Notes on tags used in MITgcmUV
     ==============================
 
+o Add moist-static energy mode (MSE_LAYERS) to package layers.
 o migrate from CVS to Git
 o pkg/cheapaml: major changes to relaxation coefficient setting:
   - rename parameters cheapaml_taurelax & cheapaml_taurelaxocean (in days)

pkg/layers/LAYERS_OPTIONS.h

@@ -26,5 +26,8 @@ C The MNC stuff is too complicated
 C Allow use of potential density as a layering field.
 #define LAYERS_PRHO_REF
 
+C Allow use of Moist Static Energy as a coordinate (relevant in the atmosphere)
+#undef MSE_LAYERS
+
 #endif /* ALLOW_LAYERS */
 #endif /* LAYERS_OPTIONS_H */

pkg/layers/layers_calc.F

@@ -29,7 +29,9 @@ SUBROUTINE LAYERS_CALC(
 #ifdef ALLOW_GMREDI
 # include "GMREDI.h"
 #endif
-
+#ifdef LAYERS_PRHO_REF
+# include "EOS.h"
+#endif
 C !INPUT PARAMETERS:
 C     myTime :: Current time in simulation
 C     myIter :: Current iteration number
@@ -71,9 +73,9 @@ SUBROUTINE LAYERS_CALC(
       _RL layers_V  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
 # endif /* LAYERS_THICKNESS */
 #endif /* LAYERS_VFLUX */
-#ifdef LAYERS_PRHO_REF
+#if (defined LAYERS_PRHO_REF) || (defined MSE_LAYERS)
       _RL prho(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
-      _RL rhoShift
+      _RL rhoShift, conv_theta2T
       INTEGER i, j, k
 #endif
 C --  other local variables:
@@ -174,12 +176,21 @@ SUBROUTINE LAYERS_CALC(
            ENDDO
           ENDDO
           DO k = 1,Nr
-           CALL FIND_RHO_2D( 1-OLx, sNx+OLx, 1-OLy, sNy+OLy,
-     &                       layers_krho(iLa),
-     &                       theta(1-OLx,1-OLy,k,bi,bj),
-     &                       salt(1-OLx,1-OLy,k,bi,bj),
-     &                       prho(1-OLx,1-OLy,k,bi,bj),
-     &                       k, bi, bj, myThid )
+           IF (equationOfState.EQ.'LINEAR') THEN
+            CALL FIND_RHO_2D( 1-OLx, sNx+OLx, 1-OLy, sNy+OLy,
+     &                        k,
+     &                        theta(1-OLx,1-OLy,k,bi,bj),
+     &                        salt(1-OLx,1-OLy,k,bi,bj),
+     &                        prho(1-OLx,1-OLy,k,bi,bj),
+     &                        k, bi, bj, myThid )
+           ELSE
+            CALL FIND_RHO_2D( 1-OLx, sNx+OLx, 1-OLy, sNy+OLy,
+     &                        layers_krho(iLa),
+     &                        theta(1-OLx,1-OLy,k,bi,bj),
+     &                        salt(1-OLx,1-OLy,k,bi,bj),
+     &                        prho(1-OLx,1-OLy,k,bi,bj),
+     &                        k, bi, bj, myThid )
+           ENDIF
 #ifdef LAYERS_THERMODYNAMICS
 C -- it might be more memory efficient not to store alpha and beta
 C    but to multiply the fluxes in place here
@@ -218,6 +229,30 @@ SUBROUTINE LAYERS_CALC(
      &              myThid)
 #endif
 #endif /* LAYERS_PRHO_REF */
+       ELSEIF ( layers_num(iLa) .EQ. 4 ) THEN
+#ifdef MSE_LAYERS
+C     For layers_num(iLa) = 4, calculate the absolute temperature referenced to 1000 mb
+        DO bj=myByLo(myThid),myByHi(myThid)
+         DO bi=myBxLo(myThid),myBxHi(myThid)
+          DO k = 1,Nr
+           conv_theta2T = (rC(k)/atm_Po)**atm_kappa
+           DO j = 1-OLy,sNy+OLy
+            DO i = 1-OLx,sNx+OLx
+             prho(i,j,k,bi,bj) = atm_Cp*conv_theta2T*theta(i,j,k,bi,bj)
+     &                         +  totPhiHyd(i,j,k,bi,bj) +  phiRef(2*k)
+     &                         +  2501. _d 0 * salt(i,j,k,bi,bj)
+            ENDDO
+           ENDDO
+          ENDDO
+         ENDDO
+        ENDDO
+        CALL LAYERS_FLUXCALC( uVel,vVel, prho, iLa,
+     &              layers_UH, layers_VH,
+     &              layers_Hw, layers_Hs,
+     &              layers_PIw,layers_PIs,
+     &              layers_U,  layers_V,
+     &              myThid )
+#endif /* MSE_LAYERS */
        ENDIF
 
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

pkg/layers/layers_readparms.F

@@ -147,6 +147,7 @@ SUBROUTINE LAYERS_READPARMS( myThid )
         IF ( layers_name(iLa).EQ.'TH ' ) layers_num(iLa) = 1
         IF ( layers_name(iLa).EQ.'SLT' ) layers_num(iLa) = 2
         IF ( layers_name(iLa).EQ.'RHO' ) layers_num(iLa) = 3
+        IF ( layers_name(iLa).EQ.'MSE' ) layers_num(iLa) = 4
         IF ( layers_name(iLa).NE.' ' .AND.
      &       layers_num(iLa).EQ.0 ) THEN
           WRITE(msgBuf,'(2A,I2,3A)') 'LAYERS_READPARMS: ',
@@ -172,6 +173,16 @@ SUBROUTINE LAYERS_READPARMS( myThid )
        ENDIF
       ENDDO
 
+#ifndef LAYERS_PRHO_REF
+      IF ( LAYER_nb .EQ. 3 ) THEN
+          WRITE(msgBuf,'(2A,I4)')
+     &      'S/R LAYERS_READPARMS: ',
+     &      'Layering in density requires to define LAYERS_PRHO_REF'
+          CALL PRINT_ERROR( msgBuf, myThid )
+          STOP 'ABNORMAL END: S/R LAYERS_READPARMS'
+      ENDIF
+#endif
+
 C--   Make sure that we locally honor the global MNC on/off flag
       layers_MNC = layers_MNC .AND. useMNC
 #ifndef ALLOW_MNC