/[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.50 by jmc, Sun Jun 14 21:45:12 2009 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 -  Number of this instance of INI_FORCING
32        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
33    
34    C     !LOCAL VARIABLES:
35  C     == Local variables ==  C     == Local variables ==
 C     iG, jG - Global coordinate index  
36  C     bi,bj  - Loop counters  C     bi,bj  - Loop counters
37  C     I,J,K  C     i, j
 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         ENDDO         ENDDO
89        ENDDO        ENDDO
90    
91        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
92         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
93          DO j=1,sNy          DO j=1-Oly,sNy+Oly
94           DO i=1,sNx           DO i=1-Olx,sNx+Olx
95            distY = yC(i,j,bi,bj)/lY            IF ( doThetaClimRelax .AND.
96  C         tauX  = -tauMax*cos(2. _d 0*PI*distY)       &         ABS(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
97            tauX  = tauMax*sin(PI*distY)             lambdaThetaClimRelax(i,j,bi,bj) = 1. _d 0/tauThetaClimRelax
98            fu(i,j,bi,bj) = tauX/(delZ(1)*rhonil)            ELSE
99  C         fu(i,j,bi,bj) = -tauMax             lambdaThetaClimRelax(i,j,bi,bj) = 0. _d 0
100              ENDIF
101              IF ( doSaltClimRelax .AND.
102         &         ABS(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
103               lambdaSaltClimRelax(i,j,bi,bj) = 1. _d 0/tauSaltClimRelax
104              ELSE
105               lambdaSaltClimRelax(i,j,bi,bj) = 0. _d 0
106              ENDIF
107           ENDDO           ENDDO
108          ENDDO          ENDDO
109         ENDDO         ENDDO
110        ENDDO        ENDDO
111  C  
112        _EXCH_XY_R4(fu  , myThid )  C-    every-one waits before master thread loads from file
113        _EXCH_XY_R4(fv  , myThid )  C     this is done within IO routines => no longer needed
114    c     _BARRIER
115  CcnhDebugStarts  
116  Cdbg  WRITE(0,*) ' distY = ', distY        IF ( zonalWindFile .NE. ' '  ) THEN
117  Cdbg  WRITE(0,*) '    ly = ', lY         CALL READ_FLD_XY_RS( zonalWindFile, ' ', fu, 0, myThid )
118  Cdbg  WRITE(0,*) ' tauMax= ', tauMax        ENDIF
119  Cdbg  CALL PLOT_FIELD_XYR8( fu, 'INI_FORCING FU',1,myThid)        IF ( meridWindFile .NE. ' '  ) THEN
120  Cdbg  STOP 'INI_FORCING'         CALL READ_FLD_XY_RS( meridWindFile, ' ', fv, 0, myThid )
121  CcnhDebugEnds        ENDIF
122          IF ( surfQFile .NE. ' '  ) THEN
123           CALL READ_FLD_XY_RS( surfQFile, ' ', Qnet, 0, myThid )
124          ELSEIF ( surfQnetFile .NE. ' '  ) THEN
125           CALL READ_FLD_XY_RS( surfQnetFile, ' ', Qnet, 0, myThid )
126          ENDIF
127          IF ( EmPmRfile .NE. ' '  ) THEN
128           CALL READ_FLD_XY_RS( EmPmRfile, ' ', EmPmR, 0, myThid )
129    c      IF ( convertEmP2rUnit.EQ.mass2rUnit ) THEN
130    C-     EmPmR is now (after c59h) expressed in kg/m2/s (fresh water mass flux)
131            DO bj = myByLo(myThid), myByHi(myThid)
132             DO bi = myBxLo(myThid), myBxHi(myThid)
133              DO j=1-Oly,sNy+Oly
134               DO i=1-Olx,sNx+Olx
135                EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj)*rhoConstFresh
136               ENDDO
137              ENDDO
138             ENDDO
139            ENDDO
140    c      ENDIF
141          ENDIF
142          IF ( saltFluxFile .NE. ' '  ) THEN
143           CALL READ_FLD_XY_RS( saltFluxFile, ' ', saltFlux, 0, myThid )
144          ENDIF
145          IF ( thetaClimFile .NE. ' '  ) THEN
146           CALL READ_FLD_XY_RS( thetaClimFile, ' ', SST, 0, myThid )
147          ENDIF
148          IF ( saltClimFile .NE. ' '  ) THEN
149           CALL READ_FLD_XY_RS( saltClimFile, ' ', SSS, 0, myThid )
150          ENDIF
151          IF ( lambdaThetaFile .NE. ' '  ) THEN
152           CALL READ_FLD_XY_RS( lambdaThetaFile, ' ',
153         &  lambdaThetaClimRelax, 0, myThid )
154          ENDIF
155          IF ( lambdaSaltFile .NE. ' '  ) THEN
156           CALL READ_FLD_XY_RS( lambdaSaltFile, ' ',
157         &  lambdaSaltClimRelax, 0, myThid )
158          ENDIF
159    #ifdef SHORTWAVE_HEATING
160          IF ( surfQswFile .NE. ' '  ) THEN
161           CALL READ_FLD_XY_RS( surfQswFile, ' ', Qsw, 0, myThid )
162           IF ( surfQFile .NE. ' '  ) THEN
163    C-     Qnet is now (after c54) the net Heat Flux (including SW)
164            DO bj = myByLo(myThid), myByHi(myThid)
165             DO bi = myBxLo(myThid), myBxHi(myThid)
166              DO j=1-OLy,sNy+OLy
167               DO i=1-OLx,sNx+OLx
168                Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) + Qsw(i,j,bi,bj)
169               ENDDO
170              ENDDO
171             ENDDO
172            ENDDO
173           ENDIF
174          ENDIF
175    #endif
176    #ifdef ATMOSPHERIC_LOADING
177          IF ( pLoadFile .NE. ' '  ) THEN
178           CALL READ_FLD_XY_RS( pLoadFile, ' ', pLoad, 0, myThid )
179          ENDIF
180    #endif
181    
182          CALL EXCH_UV_XY_RS( fu,fv, .TRUE., myThid )
183          CALL EXCH_XY_RS( Qnet , myThid )
184          CALL EXCH_XY_RS( EmPmR, myThid )
185          CALL EXCH_XY_RS( saltFlux, myThid )
186          CALL EXCH_XY_RS( SST  , myThid )
187          CALL EXCH_XY_RS( SSS  , myThid )
188          CALL EXCH_XY_RS( lambdaThetaClimRelax, myThid )
189          CALL EXCH_XY_RS( lambdaSaltClimRelax , myThid )
190    #ifdef SHORTWAVE_HEATING
191          CALL EXCH_XY_RS(Qsw  , myThid )
192    #endif
193    #ifdef ATMOSPHERIC_LOADING
194          CALL EXCH_XY_RS(pLoad  , myThid )
195    C     CALL PLOT_FIELD_XYRS( pLoad, 'S/R INI_FORCING pLoad',1,myThid)
196    #endif
197    C     CALL PLOT_FIELD_XYRS( fu, 'S/R INI_FORCING FU',1,myThid)
198    C     CALL PLOT_FIELD_XYRS( fv, 'S/R INI_FORCING FV',1,myThid)
199    
200    #ifdef ATMOSPHERIC_LOADING
201          IF ( pLoadFile .NE. ' ' .AND. usingPCoords ) THEN
202    C-- This is a hack used to read phi0surf from a file (pLoadFile)
203    C          instead of computing it from bathymetry & density ref. profile.
204    C-  Ocean: The true atmospheric P-loading is not yet implemented for P-coord
205    C          (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS).
206    C-  Atmos: sometime usefull to overwrite phi0surf with fixed-in-time field
207    C          read from file (and anyway, pressure loading is meaningless here)
208            DO bj = myByLo(myThid), myByHi(myThid)
209             DO bi = myBxLo(myThid), myBxHi(myThid)
210              DO j=1-OLy,sNy+OLy
211               DO i=1-OLx,sNx+OLx
212                 phi0surf(i,j,bi,bj) = pLoad(i,j,bi,bj)
213               ENDDO
214              ENDDO
215             ENDDO
216            ENDDO
217          ENDIF
218    #endif /* ATMOSPHERIC_LOADING */
219    
220        RETURN        RETURN
221        END        END
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.50
  ViewVC Help
Powered by ViewVC 1.1.22