model/src/ini_forcing.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_forcing.F,v 1.53 2014/08/07 18:43:33 heimbach Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

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

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE INI_FORCING
C     | o Set model initial forcing fields.
C     *==========================================================*
C     \ev

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

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

C     !LOCAL VARIABLES:
C     == Local variables ==
C     bi,bj  :: Tile indices
C     i, j   :: Loop counters
      INTEGER bi, bj
      INTEGER  i, j
CEOP

C-    Initialise all arrays in common blocks
C     <-- moved to new S/R INI_FFIELDS

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          IF ( doThetaClimRelax .AND.
     &         ABS(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
           lambdaThetaClimRelax(i,j,bi,bj) = 1. _d 0/tauThetaClimRelax
          ELSE
           lambdaThetaClimRelax(i,j,bi,bj) = 0. _d 0
          ENDIF
          IF ( doSaltClimRelax .AND.
     &         ABS(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
           lambdaSaltClimRelax(i,j,bi,bj) = 1. _d 0/tauSaltClimRelax
          ELSE
           lambdaSaltClimRelax(i,j,bi,bj) = 0. _d 0
          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDDO

C-    every-one waits before master thread loads from file
C     this is done within IO routines => no longer needed
c     _BARRIER

      IF ( zonalWindFile .NE. ' '  ) THEN
       CALL READ_FLD_XY_RS( zonalWindFile, ' ', fu, 0, myThid )
      ENDIF
      IF ( meridWindFile .NE. ' '  ) THEN
       CALL READ_FLD_XY_RS( meridWindFile, ' ', fv, 0, myThid )
      ENDIF
      IF ( surfQFile .NE. ' '  ) THEN
       CALL READ_FLD_XY_RS( surfQFile, ' ', Qnet, 0, myThid )
      ELSEIF ( surfQnetFile .NE. ' '  ) THEN
       CALL READ_FLD_XY_RS( surfQnetFile, ' ', Qnet, 0, myThid )
      ENDIF
      IF ( EmPmRfile .NE. ' '  ) THEN
       CALL READ_FLD_XY_RS( EmPmRfile, ' ', EmPmR, 0, myThid )
c      IF ( convertEmP2rUnit.EQ.mass2rUnit ) THEN
C-     EmPmR is now (after c59h) expressed in kg/m2/s (fresh water mass flux)
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj)*rhoConstFresh
           ENDDO
          ENDDO
         ENDDO
        ENDDO
c      ENDIF
      ENDIF
      IF ( saltFluxFile .NE. ' '  ) THEN
       CALL READ_FLD_XY_RS( saltFluxFile, ' ', saltFlux, 0, myThid )
      ENDIF
      IF ( thetaClimFile .NE. ' '  ) THEN
       CALL READ_FLD_XY_RS( thetaClimFile, ' ', SST, 0, myThid )
      ENDIF
      IF ( saltClimFile .NE. ' '  ) THEN
       CALL READ_FLD_XY_RS( saltClimFile, ' ', SSS, 0, myThid )
      ENDIF
      IF ( lambdaThetaFile .NE. ' '  ) THEN
       CALL READ_FLD_XY_RS( lambdaThetaFile, ' ',
     &  lambdaThetaClimRelax, 0, myThid )
      ENDIF
      IF ( lambdaSaltFile .NE. ' '  ) THEN
       CALL READ_FLD_XY_RS( lambdaSaltFile, ' ',
     &  lambdaSaltClimRelax, 0, myThid )
      ENDIF
#ifdef SHORTWAVE_HEATING
      IF ( surfQswFile .NE. ' '  ) THEN
       CALL READ_FLD_XY_RS( surfQswFile, ' ', Qsw, 0, myThid )
       IF ( surfQFile .NE. ' '  ) THEN
C-     Qnet is now (after c54) the net Heat Flux (including SW)
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) + Qsw(i,j,bi,bj)
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDIF
      ENDIF
#endif
#ifdef ATMOSPHERIC_LOADING
      IF ( pLoadFile .NE. ' '  ) THEN
       CALL READ_FLD_XY_RS( pLoadFile, ' ', pLoad, 0, myThid )
      ENDIF
#endif
#ifdef ALLOW_ADDFLUID
      IF ( addMassFile .NE. ' ' ) THEN
       CALL READ_FLD_XYZ_RL( addMassFile, ' ', addMass, 0, myThid )
       CALL EXCH_XYZ_RL( addMass, myThid )
      ENDIF
#endif /* ALLOW_ADDFLUID */
#ifdef ALLOW_GEOTHERMAL_FLUX
      IF ( geothermalFile .NE. ' ' ) THEN
       CALL READ_FLD_XY_RS( geothermalFile, ' ',
     &  geothermalFlux, 0, myThid )
       CALL EXCH_XY_RS( geothermalFlux, myThid )
# ifdef ALLOW_MONITOR
       CALL MON_PRINTSTATS_RS(
     &  1,geothermalFlux,'geothermalFlux',myThid)
# endif
      ENDIF
#endif /* ALLOW_GEOTHERMAL_FLUX */

      CALL EXCH_UV_XY_RS( fu,fv, .TRUE., myThid )
      CALL EXCH_XY_RS( Qnet , myThid )
      CALL EXCH_XY_RS( EmPmR, myThid )
      CALL EXCH_XY_RS( saltFlux, myThid )
      CALL EXCH_XY_RS( SST  , myThid )
      CALL EXCH_XY_RS( SSS  , myThid )
      CALL EXCH_XY_RS( lambdaThetaClimRelax, myThid )
      CALL EXCH_XY_RS( lambdaSaltClimRelax , myThid )
#ifdef SHORTWAVE_HEATING
      CALL EXCH_XY_RS(Qsw  , myThid )
#endif
#ifdef ATMOSPHERIC_LOADING
      CALL EXCH_XY_RS(pLoad  , myThid )
C     CALL PLOT_FIELD_XYRS( pLoad, 'S/R INI_FORCING pLoad',1,myThid)
#endif
C     CALL PLOT_FIELD_XYRS( fu, 'S/R INI_FORCING FU',1,myThid)
C     CALL PLOT_FIELD_XYRS( fv, 'S/R INI_FORCING FV',1,myThid)

#ifdef ATMOSPHERIC_LOADING
      IF ( pLoadFile .NE. ' ' .AND. usingPCoords ) THEN
C-- This is a hack used to read phi0surf from a file (pLoadFile)
C          instead of computing it from bathymetry & density ref. profile.
C-  Ocean: 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).
C-  Atmos: sometime usefull to overwrite phi0surf with fixed-in-time field
C          read from file (and anyway, pressure loading is meaningless here)
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
             phi0surf(i,j,bi,bj) = pLoad(i,j,bi,bj)
           ENDDO
          ENDDO
         ENDDO
        ENDDO
      ENDIF
#endif /* ATMOSPHERIC_LOADING */

      RETURN
      END
1	jmc	1.54	C $Header: /u/gcmpack/MITgcm/model/src/ini_forcing.F,v 1.53 2014/08/07 18:43:33 heimbach Exp $
2	cnh	1.22	C $Name: $
3	cnh	1.1
4	edhill	1.26	#include "PACKAGES_CONFIG.h"
5	cnh	1.12	#include "CPP_OPTIONS.h"
6	cnh	1.1
7	cnh	1.22	CBOP
8			C !ROUTINE: INI_FORCING
9			C !INTERFACE:
10	cnh	1.1	SUBROUTINE INI_FORCING( myThid )
11	cnh	1.22
12			C !DESCRIPTION: \bv
13			C ==========================================================
14	jmc	1.42	C \| SUBROUTINE INI_FORCING
15			C \| o Set model initial forcing fields.
16	cnh	1.22	C ==========================================================
17			C \ev
18
19			C !USES:
20	adcroft	1.14	IMPLICIT NONE
21	cnh	1.1	C === Global variables ===
22			#include "SIZE.h"
23			#include "EEPARAMS.h"
24			#include "PARAMS.h"
25			#include "GRID.h"
26	jmc	1.49	#include "SURFACE.h"
27	cnh	1.1	#include "FFIELDS.h"
28
29	cnh	1.22	C !INPUT/OUTPUT PARAMETERS:
30	cnh	1.1	C == Routine arguments ==
31	jmc	1.51	C myThid :: my Thread Id number
32	cnh	1.1	INTEGER myThid
33
34	cnh	1.22	C !LOCAL VARIABLES:
35	cnh	1.1	C == Local variables ==
36	jmc	1.51	C bi,bj :: Tile indices
37			C i, j :: Loop counters
38	cnh	1.1	INTEGER bi, bj
39	jmc	1.42	INTEGER i, j
40	cnh	1.22	CEOP
41	cnh	1.1
42	jmc	1.42	C- Initialise all arrays in common blocks
43	jmc	1.54	C <-- moved to new S/R INI_FFIELDS
44	jmc	1.42
45	heimbach	1.38	DO bj = myByLo(myThid), myByHi(myThid)
46			DO bi = myBxLo(myThid), myBxHi(myThid)
47	jmc	1.52	DO j=1-OLy,sNy+OLy
48			DO i=1-OLx,sNx+OLx
49	jmc	1.42	IF ( doThetaClimRelax .AND.
50			& ABS(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
51			lambdaThetaClimRelax(i,j,bi,bj) = 1. _d 0/tauThetaClimRelax
52	heimbach	1.38	ELSE
53	jmc	1.42	lambdaThetaClimRelax(i,j,bi,bj) = 0. _d 0
54	heimbach	1.38	ENDIF
55			IF ( doSaltClimRelax .AND.
56	jmc	1.42	& ABS(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
57			lambdaSaltClimRelax(i,j,bi,bj) = 1. _d 0/tauSaltClimRelax
58	heimbach	1.38	ELSE
59	jmc	1.42	lambdaSaltClimRelax(i,j,bi,bj) = 0. _d 0
60	heimbach	1.38	ENDIF
61			ENDDO
62			ENDDO
63			ENDDO
64			ENDDO
65	jmc	1.42
66			C- every-one waits before master thread loads from file
67	jmc	1.50	C this is done within IO routines => no longer needed
68			c _BARRIER
69	jmc	1.42
70	adcroft	1.16	IF ( zonalWindFile .NE. ' ' ) THEN
71			CALL READ_FLD_XY_RS( zonalWindFile, ' ', fu, 0, myThid )
72			ENDIF
73			IF ( meridWindFile .NE. ' ' ) THEN
74			CALL READ_FLD_XY_RS( meridWindFile, ' ', fv, 0, myThid )
75			ENDIF
76			IF ( surfQFile .NE. ' ' ) THEN
77			CALL READ_FLD_XY_RS( surfQFile, ' ', Qnet, 0, myThid )
78	jmc	1.34	ELSEIF ( surfQnetFile .NE. ' ' ) THEN
79			CALL READ_FLD_XY_RS( surfQnetFile, ' ', Qnet, 0, myThid )
80	adcroft	1.16	ENDIF
81			IF ( EmPmRfile .NE. ' ' ) THEN
82			CALL READ_FLD_XY_RS( EmPmRfile, ' ', EmPmR, 0, myThid )
83	jmc	1.47	c IF ( convertEmP2rUnit.EQ.mass2rUnit ) THEN
84	jmc	1.46	C- EmPmR is now (after c59h) expressed in kg/m2/s (fresh water mass flux)
85			DO bj = myByLo(myThid), myByHi(myThid)
86			DO bi = myBxLo(myThid), myBxHi(myThid)
87	jmc	1.52	DO j=1-OLy,sNy+OLy
88			DO i=1-OLx,sNx+OLx
89	jmc	1.46	EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj)*rhoConstFresh
90			ENDDO
91			ENDDO
92			ENDDO
93			ENDDO
94	jmc	1.47	c ENDIF
95	adcroft	1.16	ENDIF
96	jmc	1.37	IF ( saltFluxFile .NE. ' ' ) THEN
97			CALL READ_FLD_XY_RS( saltFluxFile, ' ', saltFlux, 0, myThid )
98			ENDIF
99	adcroft	1.16	IF ( thetaClimFile .NE. ' ' ) THEN
100			CALL READ_FLD_XY_RS( thetaClimFile, ' ', SST, 0, myThid )
101			ENDIF
102			IF ( saltClimFile .NE. ' ' ) THEN
103			CALL READ_FLD_XY_RS( saltClimFile, ' ', SSS, 0, myThid )
104			ENDIF
105	heimbach	1.38	IF ( lambdaThetaFile .NE. ' ' ) THEN
106	jmc	1.42	CALL READ_FLD_XY_RS( lambdaThetaFile, ' ',
107	heimbach	1.38	& lambdaThetaClimRelax, 0, myThid )
108			ENDIF
109			IF ( lambdaSaltFile .NE. ' ' ) THEN
110	jmc	1.42	CALL READ_FLD_XY_RS( lambdaSaltFile, ' ',
111	heimbach	1.38	& lambdaSaltClimRelax, 0, myThid )
112			ENDIF
113	heimbach	1.18	#ifdef SHORTWAVE_HEATING
114			IF ( surfQswFile .NE. ' ' ) THEN
115			CALL READ_FLD_XY_RS( surfQswFile, ' ', Qsw, 0, myThid )
116	jmc	1.34	IF ( surfQFile .NE. ' ' ) THEN
117			C- Qnet is now (after c54) the net Heat Flux (including SW)
118	jmc	1.42	DO bj = myByLo(myThid), myByHi(myThid)
119			DO bi = myBxLo(myThid), myBxHi(myThid)
120	jmc	1.36	DO j=1-OLy,sNy+OLy
121			DO i=1-OLx,sNx+OLx
122	jmc	1.34	Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) + Qsw(i,j,bi,bj)
123	jmc	1.36	ENDDO
124	jmc	1.34	ENDDO
125			ENDDO
126	jmc	1.36	ENDDO
127	jmc	1.34	ENDIF
128	heimbach	1.18	ENDIF
129			#endif
130	mlosch	1.23	#ifdef ATMOSPHERIC_LOADING
131			IF ( pLoadFile .NE. ' ' ) THEN
132	jmc	1.44	CALL READ_FLD_XY_RS( pLoadFile, ' ', pLoad, 0, myThid )
133	mlosch	1.23	ENDIF
134			#endif
135	jmc	1.52	#ifdef ALLOW_ADDFLUID
136			IF ( addMassFile .NE. ' ' ) THEN
137			CALL READ_FLD_XYZ_RL( addMassFile, ' ', addMass, 0, myThid )
138			CALL EXCH_XYZ_RL( addMass, myThid )
139			ENDIF
140			#endif /* ALLOW_ADDFLUID */
141	heimbach	1.53	#ifdef ALLOW_GEOTHERMAL_FLUX
142			IF ( geothermalFile .NE. ' ' ) THEN
143	jmc	1.54	CALL READ_FLD_XY_RS( geothermalFile, ' ',
144	heimbach	1.53	& geothermalFlux, 0, myThid )
145			CALL EXCH_XY_RS( geothermalFlux, myThid )
146			# ifdef ALLOW_MONITOR
147			CALL MON_PRINTSTATS_RS(
148			& 1,geothermalFlux,'geothermalFlux',myThid)
149			# endif
150			ENDIF
151			#endif /* ALLOW_GEOTHERMAL_FLUX */
152	heimbach	1.40
153	jmc	1.43	CALL EXCH_UV_XY_RS( fu,fv, .TRUE., myThid )
154			CALL EXCH_XY_RS( Qnet , myThid )
155			CALL EXCH_XY_RS( EmPmR, myThid )
156			CALL EXCH_XY_RS( saltFlux, myThid )
157			CALL EXCH_XY_RS( SST , myThid )
158			CALL EXCH_XY_RS( SSS , myThid )
159			CALL EXCH_XY_RS( lambdaThetaClimRelax, myThid )
160			CALL EXCH_XY_RS( lambdaSaltClimRelax , myThid )
161	mlosch	1.23	#ifdef SHORTWAVE_HEATING
162	jmc	1.43	CALL EXCH_XY_RS(Qsw , myThid )
163	mlosch	1.23	#endif
164			#ifdef ATMOSPHERIC_LOADING
165	jmc	1.44	CALL EXCH_XY_RS(pLoad , myThid )
166			C CALL PLOT_FIELD_XYRS( pLoad, 'S/R INI_FORCING pLoad',1,myThid)
167	mlosch	1.23	#endif
168	adcroft	1.15	C CALL PLOT_FIELD_XYRS( fu, 'S/R INI_FORCING FU',1,myThid)
169			C CALL PLOT_FIELD_XYRS( fv, 'S/R INI_FORCING FV',1,myThid)
170	cnh	1.1
171	jmc	1.49	#ifdef ATMOSPHERIC_LOADING
172			IF ( pLoadFile .NE. ' ' .AND. usingPCoords ) THEN
173			C-- This is a hack used to read phi0surf from a file (pLoadFile)
174			C instead of computing it from bathymetry & density ref. profile.
175			C- Ocean: The true atmospheric P-loading is not yet implemented for P-coord
176			C (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS).
177			C- Atmos: sometime usefull to overwrite phi0surf with fixed-in-time field
178			C read from file (and anyway, pressure loading is meaningless here)
179			DO bj = myByLo(myThid), myByHi(myThid)
180			DO bi = myBxLo(myThid), myBxHi(myThid)
181			DO j=1-OLy,sNy+OLy
182			DO i=1-OLx,sNx+OLx
183			phi0surf(i,j,bi,bj) = pLoad(i,j,bi,bj)
184			ENDDO
185			ENDDO
186			ENDDO
187			ENDDO
188			ENDIF
189			#endif /* ATMOSPHERIC_LOADING */
190
191	cnh	1.1	RETURN
192			END