/[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.1.1.1 by cnh, Wed Apr 22 19:15:30 1998 UTC revision 1.51 by jmc, Thu Apr 14 21:11:21 2011 UTC
# Line 1  Line 1 
1  C $Id$  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  CEOP
41        _RL lY  
42        _RL tauX, tauMax  C-    Initialise all arrays in common blocks
   
 C--   Initialise surface bc arrays  
 C     In cartesian yc, delY and ly are meters.  
 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  
43        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
44         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
45          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
46           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
47            fu(i,j,bi,bj) = 0. _d 0            fu              (i,j,bi,bj) = 0. _d 0
48            fv(i,j,bi,bj) = 0. _d 0            fv              (i,j,bi,bj) = 0. _d 0
49              Qnet            (i,j,bi,bj) = 0. _d 0
50              EmPmR           (i,j,bi,bj) = 0. _d 0
51              saltFlux        (i,j,bi,bj) = 0. _d 0
52              SST             (i,j,bi,bj) = 0. _d 0
53              SSS             (i,j,bi,bj) = 0. _d 0
54              Qsw             (i,j,bi,bj) = 0. _d 0
55              pLoad           (i,j,bi,bj) = 0. _d 0
56              sIceLoad        (i,j,bi,bj) = 0. _d 0
57              surfaceForcingU (i,j,bi,bj) = 0. _d 0
58              surfaceForcingV (i,j,bi,bj) = 0. _d 0
59              surfaceForcingT (i,j,bi,bj) = 0. _d 0
60              surfaceForcingS (i,j,bi,bj) = 0. _d 0
61              surfaceForcingTice(i,j,bi,bj) = 0. _d 0
62    #ifndef EXCLUDE_FFIELDS_LOAD
63              taux0           (i,j,bi,bj) = 0. _d 0
64              taux1           (i,j,bi,bj) = 0. _d 0
65              tauy0           (i,j,bi,bj) = 0. _d 0
66              tauy1           (i,j,bi,bj) = 0. _d 0
67              Qnet0           (i,j,bi,bj) = 0. _d 0
68              Qnet1           (i,j,bi,bj) = 0. _d 0
69              EmPmR0          (i,j,bi,bj) = 0. _d 0
70              EmPmR1          (i,j,bi,bj) = 0. _d 0
71              saltFlux0       (i,j,bi,bj) = 0. _d 0
72              saltFlux1       (i,j,bi,bj) = 0. _d 0
73              SST0            (i,j,bi,bj) = 0. _d 0
74              SST1            (i,j,bi,bj) = 0. _d 0
75              SSS0            (i,j,bi,bj) = 0. _d 0
76              SSS1            (i,j,bi,bj) = 0. _d 0
77    #ifdef SHORTWAVE_HEATING
78              Qsw0            (i,j,bi,bj) = 0. _d 0
79              Qsw1            (i,j,bi,bj) = 0. _d 0
80    #endif
81    #ifdef ATMOSPHERIC_LOADING
82              pLoad0          (i,j,bi,bj) = 0. _d 0
83              pLoad1          (i,j,bi,bj) = 0. _d 0
84    #endif
85    #endif /* EXCLUDE_FFIELDS_LOAD */
86           ENDDO           ENDDO
87          ENDDO          ENDDO
88    #ifndef EXCLUDE_FFIELDS_LOAD
89            loadedRec(bi,bj) = 0
90    #endif /* EXCLUDE_FFIELDS_LOAD */
91         ENDDO         ENDDO
92        ENDDO        ENDDO
93    
94        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
95         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
96          DO j=1,sNy          DO j=1-Oly,sNy+Oly
97           DO i=1,sNx           DO i=1-Olx,sNx+Olx
98            distY = yC(i,j,bi,bj)/lY            IF ( doThetaClimRelax .AND.
99  C         tauX  = -tauMax*cos(2. _d 0*PI*distY)       &         ABS(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
100            tauX  = tauMax*sin(PI*distY)             lambdaThetaClimRelax(i,j,bi,bj) = 1. _d 0/tauThetaClimRelax
101            fu(i,j,bi,bj) = tauX/(delZ(1)*rhonil)            ELSE
102  C         fu(i,j,bi,bj) = -tauMax             lambdaThetaClimRelax(i,j,bi,bj) = 0. _d 0
103              ENDIF
104              IF ( doSaltClimRelax .AND.
105         &         ABS(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
106               lambdaSaltClimRelax(i,j,bi,bj) = 1. _d 0/tauSaltClimRelax
107              ELSE
108               lambdaSaltClimRelax(i,j,bi,bj) = 0. _d 0
109              ENDIF
110           ENDDO           ENDDO
111          ENDDO          ENDDO
112         ENDDO         ENDDO
113        ENDDO        ENDDO
114  C  
115        _EXCH_XY_R4(fu  , myThid )  C-    every-one waits before master thread loads from file
116        _EXCH_XY_R4(fv  , myThid )  C     this is done within IO routines => no longer needed
117    c     _BARRIER
118  CcnhDebugStarts  
119  Cdbg  WRITE(0,*) ' distY = ', distY        IF ( zonalWindFile .NE. ' '  ) THEN
120  Cdbg  WRITE(0,*) '    ly = ', lY         CALL READ_FLD_XY_RS( zonalWindFile, ' ', fu, 0, myThid )
121  Cdbg  WRITE(0,*) ' tauMax= ', tauMax        ENDIF
122  Cdbg  CALL PLOT_FIELD_XYR8( fu, 'INI_FORCING FU',1,myThid)        IF ( meridWindFile .NE. ' '  ) THEN
123  Cdbg  STOP 'INI_FORCING'         CALL READ_FLD_XY_RS( meridWindFile, ' ', fv, 0, myThid )
124  CcnhDebugEnds        ENDIF
125          IF ( surfQFile .NE. ' '  ) THEN
126           CALL READ_FLD_XY_RS( surfQFile, ' ', Qnet, 0, myThid )
127          ELSEIF ( surfQnetFile .NE. ' '  ) THEN
128           CALL READ_FLD_XY_RS( surfQnetFile, ' ', Qnet, 0, myThid )
129          ENDIF
130          IF ( EmPmRfile .NE. ' '  ) THEN
131           CALL READ_FLD_XY_RS( EmPmRfile, ' ', EmPmR, 0, myThid )
132    c      IF ( convertEmP2rUnit.EQ.mass2rUnit ) THEN
133    C-     EmPmR is now (after c59h) expressed in kg/m2/s (fresh water mass flux)
134            DO bj = myByLo(myThid), myByHi(myThid)
135             DO bi = myBxLo(myThid), myBxHi(myThid)
136              DO j=1-Oly,sNy+Oly
137               DO i=1-Olx,sNx+Olx
138                EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj)*rhoConstFresh
139               ENDDO
140              ENDDO
141             ENDDO
142            ENDDO
143    c      ENDIF
144          ENDIF
145          IF ( saltFluxFile .NE. ' '  ) THEN
146           CALL READ_FLD_XY_RS( saltFluxFile, ' ', saltFlux, 0, myThid )
147          ENDIF
148          IF ( thetaClimFile .NE. ' '  ) THEN
149           CALL READ_FLD_XY_RS( thetaClimFile, ' ', SST, 0, myThid )
150          ENDIF
151          IF ( saltClimFile .NE. ' '  ) THEN
152           CALL READ_FLD_XY_RS( saltClimFile, ' ', SSS, 0, myThid )
153          ENDIF
154          IF ( lambdaThetaFile .NE. ' '  ) THEN
155           CALL READ_FLD_XY_RS( lambdaThetaFile, ' ',
156         &  lambdaThetaClimRelax, 0, myThid )
157          ENDIF
158          IF ( lambdaSaltFile .NE. ' '  ) THEN
159           CALL READ_FLD_XY_RS( lambdaSaltFile, ' ',
160         &  lambdaSaltClimRelax, 0, myThid )
161          ENDIF
162    #ifdef SHORTWAVE_HEATING
163          IF ( surfQswFile .NE. ' '  ) THEN
164           CALL READ_FLD_XY_RS( surfQswFile, ' ', Qsw, 0, myThid )
165           IF ( surfQFile .NE. ' '  ) THEN
166    C-     Qnet is now (after c54) the net Heat Flux (including SW)
167            DO bj = myByLo(myThid), myByHi(myThid)
168             DO bi = myBxLo(myThid), myBxHi(myThid)
169              DO j=1-OLy,sNy+OLy
170               DO i=1-OLx,sNx+OLx
171                Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) + Qsw(i,j,bi,bj)
172               ENDDO
173              ENDDO
174             ENDDO
175            ENDDO
176           ENDIF
177          ENDIF
178    #endif
179    #ifdef ATMOSPHERIC_LOADING
180          IF ( pLoadFile .NE. ' '  ) THEN
181           CALL READ_FLD_XY_RS( pLoadFile, ' ', pLoad, 0, myThid )
182          ENDIF
183    #endif
184    
185          CALL EXCH_UV_XY_RS( fu,fv, .TRUE., myThid )
186          CALL EXCH_XY_RS( Qnet , myThid )
187          CALL EXCH_XY_RS( EmPmR, myThid )
188          CALL EXCH_XY_RS( saltFlux, myThid )
189          CALL EXCH_XY_RS( SST  , myThid )
190          CALL EXCH_XY_RS( SSS  , myThid )
191          CALL EXCH_XY_RS( lambdaThetaClimRelax, myThid )
192          CALL EXCH_XY_RS( lambdaSaltClimRelax , myThid )
193    #ifdef SHORTWAVE_HEATING
194          CALL EXCH_XY_RS(Qsw  , myThid )
195    #endif
196    #ifdef ATMOSPHERIC_LOADING
197          CALL EXCH_XY_RS(pLoad  , myThid )
198    C     CALL PLOT_FIELD_XYRS( pLoad, 'S/R INI_FORCING pLoad',1,myThid)
199    #endif
200    C     CALL PLOT_FIELD_XYRS( fu, 'S/R INI_FORCING FU',1,myThid)
201    C     CALL PLOT_FIELD_XYRS( fv, 'S/R INI_FORCING FV',1,myThid)
202    
203    #ifdef ATMOSPHERIC_LOADING
204          IF ( pLoadFile .NE. ' ' .AND. usingPCoords ) THEN
205    C-- This is a hack used to read phi0surf from a file (pLoadFile)
206    C          instead of computing it from bathymetry & density ref. profile.
207    C-  Ocean: The true atmospheric P-loading is not yet implemented for P-coord
208    C          (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS).
209    C-  Atmos: sometime usefull to overwrite phi0surf with fixed-in-time field
210    C          read from file (and anyway, pressure loading is meaningless here)
211            DO bj = myByLo(myThid), myByHi(myThid)
212             DO bi = myBxLo(myThid), myBxHi(myThid)
213              DO j=1-OLy,sNy+OLy
214               DO i=1-OLx,sNx+OLx
215                 phi0surf(i,j,bi,bj) = pLoad(i,j,bi,bj)
216               ENDDO
217              ENDDO
218             ENDDO
219            ENDDO
220          ENDIF
221    #endif /* ATMOSPHERIC_LOADING */
222    
223        RETURN        RETURN
224        END        END
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.51
  ViewVC Help
Powered by ViewVC 1.1.22