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	C $Header: /u/gcmpack/MITgcm/model/src/ini_forcing.F,v 1.53 2014/08/07 18:43:33 heimbach Exp $
2	C $Name: $
3
4	#include "PACKAGES_CONFIG.h"
5	#include "CPP_OPTIONS.h"
6
7	CBOP
8	C !ROUTINE: INI_FORCING
9	C !INTERFACE:
10	SUBROUTINE INI_FORCING( myThid )
11
12	C !DESCRIPTION: \bv
13	C ==========================================================
14	C \| SUBROUTINE INI_FORCING
15	C \| o Set model initial forcing fields.
16	C ==========================================================
17	C \ev
18
19	C !USES:
20	IMPLICIT NONE
21	C === Global variables ===
22	#include "SIZE.h"
23	#include "EEPARAMS.h"
24	#include "PARAMS.h"
25	#include "GRID.h"
26	#include "SURFACE.h"
27	#include "FFIELDS.h"
28
29	C !INPUT/OUTPUT PARAMETERS:
30	C == Routine arguments ==
31	C myThid :: my Thread Id number
32	INTEGER myThid
33
34	C !LOCAL VARIABLES:
35	C == Local variables ==
36	C bi,bj :: Tile indices
37	C i, j :: Loop counters
38	INTEGER bi, bj
39	INTEGER i, j
40	CEOP
41
42	C- Initialise all arrays in common blocks
43	C <-- moved to new S/R INI_FFIELDS
44
45	DO bj = myByLo(myThid), myByHi(myThid)
46	DO bi = myBxLo(myThid), myBxHi(myThid)
47	DO j=1-OLy,sNy+OLy
48	DO i=1-OLx,sNx+OLx
49	IF ( doThetaClimRelax .AND.
50	& ABS(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
51	lambdaThetaClimRelax(i,j,bi,bj) = 1. _d 0/tauThetaClimRelax
52	ELSE
53	lambdaThetaClimRelax(i,j,bi,bj) = 0. _d 0
54	ENDIF
55	IF ( doSaltClimRelax .AND.
56	& ABS(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
57	lambdaSaltClimRelax(i,j,bi,bj) = 1. _d 0/tauSaltClimRelax
58	ELSE
59	lambdaSaltClimRelax(i,j,bi,bj) = 0. _d 0
60	ENDIF
61	ENDDO
62	ENDDO
63	ENDDO
64	ENDDO
65
66	C- every-one waits before master thread loads from file
67	C this is done within IO routines => no longer needed
68	c _BARRIER
69
70	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	ELSEIF ( surfQnetFile .NE. ' ' ) THEN
79	CALL READ_FLD_XY_RS( surfQnetFile, ' ', Qnet, 0, myThid )
80	ENDIF
81	IF ( EmPmRfile .NE. ' ' ) THEN
82	CALL READ_FLD_XY_RS( EmPmRfile, ' ', EmPmR, 0, myThid )
83	c IF ( convertEmP2rUnit.EQ.mass2rUnit ) THEN
84	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	DO j=1-OLy,sNy+OLy
88	DO i=1-OLx,sNx+OLx
89	EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj)*rhoConstFresh
90	ENDDO
91	ENDDO
92	ENDDO
93	ENDDO
94	c ENDIF
95	ENDIF
96	IF ( saltFluxFile .NE. ' ' ) THEN
97	CALL READ_FLD_XY_RS( saltFluxFile, ' ', saltFlux, 0, myThid )
98	ENDIF
99	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	IF ( lambdaThetaFile .NE. ' ' ) THEN
106	CALL READ_FLD_XY_RS( lambdaThetaFile, ' ',
107	& lambdaThetaClimRelax, 0, myThid )
108	ENDIF
109	IF ( lambdaSaltFile .NE. ' ' ) THEN
110	CALL READ_FLD_XY_RS( lambdaSaltFile, ' ',
111	& lambdaSaltClimRelax, 0, myThid )
112	ENDIF
113	#ifdef SHORTWAVE_HEATING
114	IF ( surfQswFile .NE. ' ' ) THEN
115	CALL READ_FLD_XY_RS( surfQswFile, ' ', Qsw, 0, myThid )
116	IF ( surfQFile .NE. ' ' ) THEN
117	C- Qnet is now (after c54) the net Heat Flux (including SW)
118	DO bj = myByLo(myThid), myByHi(myThid)
119	DO bi = myBxLo(myThid), myBxHi(myThid)
120	DO j=1-OLy,sNy+OLy
121	DO i=1-OLx,sNx+OLx
122	Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) + Qsw(i,j,bi,bj)
123	ENDDO
124	ENDDO
125	ENDDO
126	ENDDO
127	ENDIF
128	ENDIF
129	#endif
130	#ifdef ATMOSPHERIC_LOADING
131	IF ( pLoadFile .NE. ' ' ) THEN
132	CALL READ_FLD_XY_RS( pLoadFile, ' ', pLoad, 0, myThid )
133	ENDIF
134	#endif
135	#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	#ifdef ALLOW_GEOTHERMAL_FLUX
142	IF ( geothermalFile .NE. ' ' ) THEN
143	CALL READ_FLD_XY_RS( geothermalFile, ' ',
144	& 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
153	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	#ifdef SHORTWAVE_HEATING
162	CALL EXCH_XY_RS(Qsw , myThid )
163	#endif
164	#ifdef ATMOSPHERIC_LOADING
165	CALL EXCH_XY_RS(pLoad , myThid )
166	C CALL PLOT_FIELD_XYRS( pLoad, 'S/R INI_FORCING pLoad',1,myThid)
167	#endif
168	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
171	#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	RETURN
192	END