/[MITgcm]/MITgcm/model/src/ini_forcing.F
ViewVC logotype

Diff of /MITgcm/model/src/ini_forcing.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.2 by cnh, Fri Apr 24 02:05:41 1998 UTC revision 1.52 by jmc, Fri Nov 9 22:43:53 2012 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "PACKAGES_CONFIG.h"
5    #include "CPP_OPTIONS.h"
6    
7  CStartOfInterface  CBOP
8    C     !ROUTINE: INI_FORCING
9    C     !INTERFACE:
10        SUBROUTINE INI_FORCING( myThid )        SUBROUTINE INI_FORCING( myThid )
 C     /==========================================================\  
 C     | SUBROUTINE INI_FORCING                                   |  
 C     | o Set model initial forcing fields.                      |  
 C     \==========================================================/  
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 ===  C     === Global variables ===
22  #include "SIZE.h"  #include "SIZE.h"
23  #include "EEPARAMS.h"  #include "EEPARAMS.h"
24  #include "PARAMS.h"  #include "PARAMS.h"
25  #include "GRID.h"  #include "GRID.h"
26    #include "SURFACE.h"
27  #include "FFIELDS.h"  #include "FFIELDS.h"
28    
29    C     !INPUT/OUTPUT PARAMETERS:
30  C     == Routine arguments ==  C     == Routine arguments ==
31  C     myThid -  Number of this instance of INI_FORCING  C     myThid :: my Thread Id number
32        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
33    
34    C     !LOCAL VARIABLES:
35  C     == Local variables ==  C     == Local variables ==
36  C     iG, jG - Global coordinate index  C     bi,bj  :: Tile indices
37  C     bi,bj  - Loop counters  C     i, j   :: Loop counters
 C     I,J,K  
 C     distY  - Distance accross domain of point j  
 C     lY     - Width of the basin ( last point is land )  
 C     tauMax - Peak zonal wind stress ( N/m^2 )  
 C     tauX   - Local zonal wind stress  
       INTEGER iG, jG  
38        INTEGER bi, bj        INTEGER bi, bj
39        INTEGER  I,  J, K        INTEGER  i, j
40        _RL distY  #if (defined ALLOW_ADDFLUID) || (defined ALLOW_FRICTION_HEATING)
41        _RL lY        INTEGER k
42        _RL tauX, tauMax  #endif
43    CEOP
44  C--   Initialise surface bc arrays  
45  C     In cartesian yc, delY and ly are meters.  C-    Initialise all arrays in common blocks
 C     In spherical polar yc, delY and ly are degrees  
       tauMax = 0.1 _d 0  
       lY = 0. _d 0  
       DO j=1,nY-1  
        lY = lY + delY(j)  
       ENDDO  
46        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
47         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
48          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
49           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
50            fu(i,j,bi,bj) = 0. _d 0            fu              (i,j,bi,bj) = 0. _d 0
51            fv(i,j,bi,bj) = 0. _d 0            fv              (i,j,bi,bj) = 0. _d 0
52              Qnet            (i,j,bi,bj) = 0. _d 0
53              EmPmR           (i,j,bi,bj) = 0. _d 0
54              saltFlux        (i,j,bi,bj) = 0. _d 0
55              SST             (i,j,bi,bj) = 0. _d 0
56              SSS             (i,j,bi,bj) = 0. _d 0
57              Qsw             (i,j,bi,bj) = 0. _d 0
58              pLoad           (i,j,bi,bj) = 0. _d 0
59              sIceLoad        (i,j,bi,bj) = 0. _d 0
60              surfaceForcingU (i,j,bi,bj) = 0. _d 0
61              surfaceForcingV (i,j,bi,bj) = 0. _d 0
62              surfaceForcingT (i,j,bi,bj) = 0. _d 0
63              surfaceForcingS (i,j,bi,bj) = 0. _d 0
64              surfaceForcingTice(i,j,bi,bj) = 0. _d 0
65    #ifndef EXCLUDE_FFIELDS_LOAD
66              taux0           (i,j,bi,bj) = 0. _d 0
67              taux1           (i,j,bi,bj) = 0. _d 0
68              tauy0           (i,j,bi,bj) = 0. _d 0
69              tauy1           (i,j,bi,bj) = 0. _d 0
70              Qnet0           (i,j,bi,bj) = 0. _d 0
71              Qnet1           (i,j,bi,bj) = 0. _d 0
72              EmPmR0          (i,j,bi,bj) = 0. _d 0
73              EmPmR1          (i,j,bi,bj) = 0. _d 0
74              saltFlux0       (i,j,bi,bj) = 0. _d 0
75              saltFlux1       (i,j,bi,bj) = 0. _d 0
76              SST0            (i,j,bi,bj) = 0. _d 0
77              SST1            (i,j,bi,bj) = 0. _d 0
78              SSS0            (i,j,bi,bj) = 0. _d 0
79              SSS1            (i,j,bi,bj) = 0. _d 0
80    #ifdef SHORTWAVE_HEATING
81              Qsw0            (i,j,bi,bj) = 0. _d 0
82              Qsw1            (i,j,bi,bj) = 0. _d 0
83    #endif
84    #ifdef ATMOSPHERIC_LOADING
85              pLoad0          (i,j,bi,bj) = 0. _d 0
86              pLoad1          (i,j,bi,bj) = 0. _d 0
87    #endif
88    #endif /* EXCLUDE_FFIELDS_LOAD */
89             ENDDO
90            ENDDO
91    #ifndef EXCLUDE_FFIELDS_LOAD
92            loadedRec(bi,bj) = 0
93    #endif /* EXCLUDE_FFIELDS_LOAD */
94            DO j=1-OLy,sNy+OLy
95             DO i=1-OLx,sNx+OLx
96             ENDDO
97            ENDDO
98    #ifdef ALLOW_ADDFLUID
99            DO k=1,Nr
100             DO j=1-OLy,sNy+OLy
101              DO i=1-OLx,sNx+OLx
102               addMass(i,j,k,bi,bj) = 0. _d 0
103              ENDDO
104             ENDDO
105            ENDDO
106    #endif /* ALLOW_ADDFLUID */
107    #ifdef ALLOW_FRICTION_HEATING
108            DO k=1,Nr
109             DO j=1-OLy,sNy+OLy
110              DO i=1-OLx,sNx+OLx
111               frictionHeating(i,j,k,bi,bj) = 0. _d 0
112              ENDDO
113           ENDDO           ENDDO
114          ENDDO          ENDDO
115    #endif /* ALLOW_FRICTION_HEATING */
116         ENDDO         ENDDO
117        ENDDO        ENDDO
118    
119        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
120         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
121          DO j=1,sNy          DO j=1-OLy,sNy+OLy
122           DO i=1,sNx           DO i=1-OLx,sNx+OLx
123            distY = yC(i,j,bi,bj)/lY            IF ( doThetaClimRelax .AND.
124  C         tauX  = -tauMax*cos(2. _d 0*PI*distY)       &         ABS(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
125            tauX  = tauMax*sin(PI*distY)             lambdaThetaClimRelax(i,j,bi,bj) = 1. _d 0/tauThetaClimRelax
126            fu(i,j,bi,bj) = tauX/(delZ(1)*rhonil)            ELSE
127  C         fu(i,j,bi,bj) = -tauMax             lambdaThetaClimRelax(i,j,bi,bj) = 0. _d 0
128              ENDIF
129              IF ( doSaltClimRelax .AND.
130         &         ABS(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
131               lambdaSaltClimRelax(i,j,bi,bj) = 1. _d 0/tauSaltClimRelax
132              ELSE
133               lambdaSaltClimRelax(i,j,bi,bj) = 0. _d 0
134              ENDIF
135           ENDDO           ENDDO
136          ENDDO          ENDDO
137         ENDDO         ENDDO
138        ENDDO        ENDDO
139  C  
140        _EXCH_XY_R4(fu  , myThid )  C-    every-one waits before master thread loads from file
141        _EXCH_XY_R4(fv  , myThid )  C     this is done within IO routines => no longer needed
142    c     _BARRIER
143  CcnhDebugStarts  
144  Cdbg  WRITE(0,*) ' distY = ', distY        IF ( zonalWindFile .NE. ' '  ) THEN
145  Cdbg  WRITE(0,*) '    ly = ', lY         CALL READ_FLD_XY_RS( zonalWindFile, ' ', fu, 0, myThid )
146  Cdbg  WRITE(0,*) ' tauMax= ', tauMax        ENDIF
147  Cdbg  CALL PLOT_FIELD_XYR8( fu, 'INI_FORCING FU',1,myThid)        IF ( meridWindFile .NE. ' '  ) THEN
148  Cdbg  STOP 'INI_FORCING'         CALL READ_FLD_XY_RS( meridWindFile, ' ', fv, 0, myThid )
149  CcnhDebugEnds        ENDIF
150          IF ( surfQFile .NE. ' '  ) THEN
151           CALL READ_FLD_XY_RS( surfQFile, ' ', Qnet, 0, myThid )
152          ELSEIF ( surfQnetFile .NE. ' '  ) THEN
153           CALL READ_FLD_XY_RS( surfQnetFile, ' ', Qnet, 0, myThid )
154          ENDIF
155          IF ( EmPmRfile .NE. ' '  ) THEN
156           CALL READ_FLD_XY_RS( EmPmRfile, ' ', EmPmR, 0, myThid )
157    c      IF ( convertEmP2rUnit.EQ.mass2rUnit ) THEN
158    C-     EmPmR is now (after c59h) expressed in kg/m2/s (fresh water mass flux)
159            DO bj = myByLo(myThid), myByHi(myThid)
160             DO bi = myBxLo(myThid), myBxHi(myThid)
161              DO j=1-OLy,sNy+OLy
162               DO i=1-OLx,sNx+OLx
163                EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj)*rhoConstFresh
164               ENDDO
165              ENDDO
166             ENDDO
167            ENDDO
168    c      ENDIF
169          ENDIF
170          IF ( saltFluxFile .NE. ' '  ) THEN
171           CALL READ_FLD_XY_RS( saltFluxFile, ' ', saltFlux, 0, myThid )
172          ENDIF
173          IF ( thetaClimFile .NE. ' '  ) THEN
174           CALL READ_FLD_XY_RS( thetaClimFile, ' ', SST, 0, myThid )
175          ENDIF
176          IF ( saltClimFile .NE. ' '  ) THEN
177           CALL READ_FLD_XY_RS( saltClimFile, ' ', SSS, 0, myThid )
178          ENDIF
179          IF ( lambdaThetaFile .NE. ' '  ) THEN
180           CALL READ_FLD_XY_RS( lambdaThetaFile, ' ',
181         &  lambdaThetaClimRelax, 0, myThid )
182          ENDIF
183          IF ( lambdaSaltFile .NE. ' '  ) THEN
184           CALL READ_FLD_XY_RS( lambdaSaltFile, ' ',
185         &  lambdaSaltClimRelax, 0, myThid )
186          ENDIF
187    #ifdef SHORTWAVE_HEATING
188          IF ( surfQswFile .NE. ' '  ) THEN
189           CALL READ_FLD_XY_RS( surfQswFile, ' ', Qsw, 0, myThid )
190           IF ( surfQFile .NE. ' '  ) THEN
191    C-     Qnet is now (after c54) the net Heat Flux (including SW)
192            DO bj = myByLo(myThid), myByHi(myThid)
193             DO bi = myBxLo(myThid), myBxHi(myThid)
194              DO j=1-OLy,sNy+OLy
195               DO i=1-OLx,sNx+OLx
196                Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) + Qsw(i,j,bi,bj)
197               ENDDO
198              ENDDO
199             ENDDO
200            ENDDO
201           ENDIF
202          ENDIF
203    #endif
204    #ifdef ATMOSPHERIC_LOADING
205          IF ( pLoadFile .NE. ' '  ) THEN
206           CALL READ_FLD_XY_RS( pLoadFile, ' ', pLoad, 0, myThid )
207          ENDIF
208    #endif
209    #ifdef ALLOW_ADDFLUID
210          IF ( addMassFile .NE. ' ' ) THEN
211           CALL READ_FLD_XYZ_RL( addMassFile, ' ', addMass, 0, myThid )
212           CALL EXCH_XYZ_RL( addMass, myThid )
213          ENDIF
214    #endif /* ALLOW_ADDFLUID */
215    
216          CALL EXCH_UV_XY_RS( fu,fv, .TRUE., myThid )
217          CALL EXCH_XY_RS( Qnet , myThid )
218          CALL EXCH_XY_RS( EmPmR, myThid )
219          CALL EXCH_XY_RS( saltFlux, myThid )
220          CALL EXCH_XY_RS( SST  , myThid )
221          CALL EXCH_XY_RS( SSS  , myThid )
222          CALL EXCH_XY_RS( lambdaThetaClimRelax, myThid )
223          CALL EXCH_XY_RS( lambdaSaltClimRelax , myThid )
224    #ifdef SHORTWAVE_HEATING
225          CALL EXCH_XY_RS(Qsw  , myThid )
226    #endif
227    #ifdef ATMOSPHERIC_LOADING
228          CALL EXCH_XY_RS(pLoad  , myThid )
229    C     CALL PLOT_FIELD_XYRS( pLoad, 'S/R INI_FORCING pLoad',1,myThid)
230    #endif
231    C     CALL PLOT_FIELD_XYRS( fu, 'S/R INI_FORCING FU',1,myThid)
232    C     CALL PLOT_FIELD_XYRS( fv, 'S/R INI_FORCING FV',1,myThid)
233    
234    #ifdef ATMOSPHERIC_LOADING
235          IF ( pLoadFile .NE. ' ' .AND. usingPCoords ) THEN
236    C-- This is a hack used to read phi0surf from a file (pLoadFile)
237    C          instead of computing it from bathymetry & density ref. profile.
238    C-  Ocean: The true atmospheric P-loading is not yet implemented for P-coord
239    C          (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS).
240    C-  Atmos: sometime usefull to overwrite phi0surf with fixed-in-time field
241    C          read from file (and anyway, pressure loading is meaningless here)
242            DO bj = myByLo(myThid), myByHi(myThid)
243             DO bi = myBxLo(myThid), myBxHi(myThid)
244              DO j=1-OLy,sNy+OLy
245               DO i=1-OLx,sNx+OLx
246                 phi0surf(i,j,bi,bj) = pLoad(i,j,bi,bj)
247               ENDDO
248              ENDDO
249             ENDDO
250            ENDDO
251          ENDIF
252    #endif /* ATMOSPHERIC_LOADING */
253    
254        RETURN        RETURN
255        END        END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.52
  ViewVC Help
Powered by ViewVC 1.1.22