model/src/ini_nh_fields.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_nh_fields.F,v 1.5 2013/07/28 21:06:00 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: INI_NH_FIELDS
C     !INTERFACE:
      SUBROUTINE INI_NH_FIELDS( myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE INI_NH_FIELDS
C     | o Set model initial non-hydrostatic fields.
C     *==========================================================*
C     | Note: If using NH form,
C     |  call this S/R whether starting or restarting simulation.
C     | This is different from other "true" ini_fields type S/R
C     |  (e.g., INI_VEL) which are called only when starting.
C     | Reason: no real physical field to initialise (since wVel
C     |  is diagnose from continuity) but needs to set few arrays
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "RESTART.h"
#include "NH_VARS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myThid  :: My Thread Id number
      INTEGER myThid

#ifdef ALLOW_NONHYDROSTATIC
C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER bi,bj
      INTEGER i,j
      INTEGER ks
      CHARACTER*(MAX_LEN_MBUF) msgBuf
#ifdef NONLIN_FRSURF
      INTEGER k
#endif
CEOP

      IF ( startTime .EQ. baseTime .AND.  nIter0 .EQ. 0
     &     .AND. pickupSuff .EQ. ' ' ) THEN
C--   Case where starting from initial conditions

C--   Read an initial non-hydrostatic pressure field
c       IF (phiNHinitFile .NE. ' ') THEN
c        CALL READ_FLD_XYZ_RL( phiNHinitFile, ' ', phi_nh, 0, myThid )
c        _EXCH_XYZ_RL(phi_nh, myThid)
c       ENDIF

      ELSE
C--   Case where restarting from a pickup

       _BEGIN_MASTER(myThid)
       WRITE(msgBuf,'(A,I4)')
     &   'INI_NH_FIELDS: dPhiNHstatus=', dPhiNHstatus
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                     SQUEEZE_RIGHT, myThid )
       _END_MASTER(myThid)
       IF ( exactConserv .AND. dPhiNHstatus.EQ.0 ) THEN
c      IF ( exactConserv ) THEN
C--   Separate the Hydrostatic Surface Pressure adjusment (=> put it in dPhiNH)
C     from the Non-hydrostatic pressure (since cg3d_x contains both contribution)
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          IF ( select_rStar.EQ.0 .AND. uniformFreeSurfLev ) THEN
C-       Z coordinate: assume surface @ level k=1
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
              dPhiNH(i,j,bi,bj) = phi_nh(i,j,1,bi,bj)
c             dPhiNH(i,j,bi,bj) = 0.
            ENDDO
           ENDDO
          ELSEIF ( select_rStar.EQ.0 ) THEN
C-       Other than Z coordinate: no assumption on surface level index
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             ks = kSurfC(i,j,bi,bj)
             IF ( ks.LE.Nr ) THEN
              dPhiNH(i,j,bi,bj) = phi_nh(i,j,ks,bi,bj)
             ELSE
              dPhiNH(i,j,bi,bj) = 0.
             ENDIF
            ENDDO
           ENDDO
#ifdef NONLIN_FRSURF
          ELSE
C        rStar : take vertical average of P_NH as Hyd.Surf.Press adjustment
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
              dPhiNH(i,j,bi,bj) = 0.
            ENDDO
           ENDDO
           DO k=1,Nr
            DO j=1-OLy,sNy+OLy
             DO i=1-OLx,sNx+OLx
              dPhiNH(i,j,bi,bj) = dPhiNH(i,j,bi,bj)
     &          + phi_nh(i,j,k,bi,bj)*drF(k)*hFacC(i,j,k,bi,bj)
             ENDDO
            ENDDO
           ENDDO
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
              dPhiNH(i,j,bi,bj) = dPhiNH(i,j,bi,bj)
     &                           *recip_Rcol(i,j,bi,bj)
            ENDDO
           ENDDO
#endif /* NONLIN_FRSURF */
          ENDIF
         ENDDO
        ENDDO
C-     end of if-block: dPhiNH_status
       ENDIF

      ENDIF

#endif /* ALLOW_NONHYDROSTATIC */
      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm/model/src/ini_nh_fields.F,v 1.5 2013/07/28 21:06:00 jmc Exp $
2	C $Name: $
3
4	#include "CPP_OPTIONS.h"
5
6	CBOP
7	C !ROUTINE: INI_NH_FIELDS
8	C !INTERFACE:
9	SUBROUTINE INI_NH_FIELDS( myThid )
10
11	C !DESCRIPTION: \bv
12	C ==========================================================
13	C \| SUBROUTINE INI_NH_FIELDS
14	C \| o Set model initial non-hydrostatic fields.
15	C ==========================================================
16	C \| Note: If using NH form,
17	C \| call this S/R whether starting or restarting simulation.
18	C \| This is different from other "true" ini_fields type S/R
19	C \| (e.g., INI_VEL) which are called only when starting.
20	C \| Reason: no real physical field to initialise (since wVel
21	C \| is diagnose from continuity) but needs to set few arrays
22	C ==========================================================
23	C \ev
24
25	C !USES:
26	IMPLICIT NONE
27	C === Global variables ===
28	#include "SIZE.h"
29	#include "EEPARAMS.h"
30	#include "PARAMS.h"
31	#include "GRID.h"
32	#include "RESTART.h"
33	#include "NH_VARS.h"
34
35	C !INPUT/OUTPUT PARAMETERS:
36	C == Routine arguments ==
37	C myThid :: My Thread Id number
38	INTEGER myThid
39
40	#ifdef ALLOW_NONHYDROSTATIC
41	C !LOCAL VARIABLES:
42	C == Local variables ==
43	INTEGER bi,bj
44	INTEGER i,j
45	INTEGER ks
46	CHARACTER*(MAX_LEN_MBUF) msgBuf
47	#ifdef NONLIN_FRSURF
48	INTEGER k
49	#endif
50	CEOP
51
52	IF ( startTime .EQ. baseTime .AND. nIter0 .EQ. 0
53	& .AND. pickupSuff .EQ. ' ' ) THEN
54	C-- Case where starting from initial conditions
55
56	C-- Read an initial non-hydrostatic pressure field
57	c IF (phiNHinitFile .NE. ' ') THEN
58	c CALL READ_FLD_XYZ_RL( phiNHinitFile, ' ', phi_nh, 0, myThid )
59	c _EXCH_XYZ_RL(phi_nh, myThid)
60	c ENDIF
61
62	ELSE
63	C-- Case where restarting from a pickup
64
65	_BEGIN_MASTER(myThid)
66	WRITE(msgBuf,'(A,I4)')
67	& 'INI_NH_FIELDS: dPhiNHstatus=', dPhiNHstatus
68	CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
69	& SQUEEZE_RIGHT, myThid )
70	_END_MASTER(myThid)
71	IF ( exactConserv .AND. dPhiNHstatus.EQ.0 ) THEN
72	c IF ( exactConserv ) THEN
73	C-- Separate the Hydrostatic Surface Pressure adjusment (=> put it in dPhiNH)
74	C from the Non-hydrostatic pressure (since cg3d_x contains both contribution)
75	DO bj=myByLo(myThid),myByHi(myThid)
76	DO bi=myBxLo(myThid),myBxHi(myThid)
77	IF ( select_rStar.EQ.0 .AND. uniformFreeSurfLev ) THEN
78	C- Z coordinate: assume surface @ level k=1
79	DO j=1-OLy,sNy+OLy
80	DO i=1-OLx,sNx+OLx
81	dPhiNH(i,j,bi,bj) = phi_nh(i,j,1,bi,bj)
82	c dPhiNH(i,j,bi,bj) = 0.
83	ENDDO
84	ENDDO
85	ELSEIF ( select_rStar.EQ.0 ) THEN
86	C- Other than Z coordinate: no assumption on surface level index
87	DO j=1-OLy,sNy+OLy
88	DO i=1-OLx,sNx+OLx
89	ks = kSurfC(i,j,bi,bj)
90	IF ( ks.LE.Nr ) THEN
91	dPhiNH(i,j,bi,bj) = phi_nh(i,j,ks,bi,bj)
92	ELSE
93	dPhiNH(i,j,bi,bj) = 0.
94	ENDIF
95	ENDDO
96	ENDDO
97	#ifdef NONLIN_FRSURF
98	ELSE
99	C rStar : take vertical average of P_NH as Hyd.Surf.Press adjustment
100	DO j=1-OLy,sNy+OLy
101	DO i=1-OLx,sNx+OLx
102	dPhiNH(i,j,bi,bj) = 0.
103	ENDDO
104	ENDDO
105	DO k=1,Nr
106	DO j=1-OLy,sNy+OLy
107	DO i=1-OLx,sNx+OLx
108	dPhiNH(i,j,bi,bj) = dPhiNH(i,j,bi,bj)
109	& + phi_nh(i,j,k,bi,bj)drF(k)hFacC(i,j,k,bi,bj)
110	ENDDO
111	ENDDO
112	ENDDO
113	DO j=1-OLy,sNy+OLy
114	DO i=1-OLx,sNx+OLx
115	dPhiNH(i,j,bi,bj) = dPhiNH(i,j,bi,bj)
116	& *recip_Rcol(i,j,bi,bj)
117	ENDDO
118	ENDDO
119	#endif /* NONLIN_FRSURF */
120	ENDIF
121	ENDDO
122	ENDDO
123	C- end of if-block: dPhiNH_status
124	ENDIF
125
126	ENDIF
127
128	#endif /* ALLOW_NONHYDROSTATIC */
129	RETURN
130	END