/[MITgcm]/MITgcm/pkg/cheapaml/cheapaml_fields_load.F
ViewVC logotype

Diff of /MITgcm/pkg/cheapaml/cheapaml_fields_load.F

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

revision 1.7 by jmc, Sun Feb 6 00:18:05 2011 UTC revision 1.8 by wienders, Thu Feb 24 16:11:41 2011 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "CHEAPAML_OPTIONS.h"  #include "CHEAPAML_OPTIONS.h"
5    
6    
7  C     !ROUTINE: CHEAPAML_FIELDS_LOAD  C     !ROUTINE: CHEAPAML_FIELDS_LOAD
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE CHEAPAML_FIELDS_LOAD( myTime, myIter, myThid )        SUBROUTINE CHEAPAML_FIELDS_LOAD( myTime, myIter, myThid )
# Line 21  C     === Global variables === Line 22  C     === Global variables ===
22  c #include "GRID.h"  c #include "GRID.h"
23  c #include "DYNVARS.h"  c #include "DYNVARS.h"
24  C #include "BULKF.h"  C #include "BULKF.h"
25  c #ifdef ALLOW_THSICE  c#ifdef ALLOW_THSICE
26  c #include "THSICE_VARS.h"  c#include "THSICE_VARS.h"
27  c #endif  c#endif
28  #include "CHEAPAML.h"  #include "CHEAPAML.h"
29    
30  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
# Line 64  C     aWght, bWght :: Interpolation weig Line 65  C     aWght, bWght :: Interpolation weig
65       &         wavesp0,  wavesp1,       &         wavesp0,  wavesp1,
66       &         rair0,    rair1       &         rair0,    rair1
67    
68    
69        _RL  trair0    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  trair0    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
70        _RL  trair1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  trair1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
71        _RL  qrair0    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  qrair0    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
# Line 85  C     aWght, bWght :: Interpolation weig Line 87  C     aWght, bWght :: Interpolation weig
87        _RL  rair0  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  rair0  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
88        _RL  rair1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  rair1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
89    
90        INTEGER bi,bj,i,j,intime0,intime1          INTEGER bi,bj,i,j,intime0,intime1
91        INTEGER iG,jG          INTEGER iG,jG
92        _RL aWght,bWght,rdt,u          _RL aWght,bWght,rdt,u
93        _RL ssq0,ssq1,ssq2,ssqa          _RL ssq0,ssq1,ssq2,ssqa
94  c xsolph - phase of year, assuming time zero is mid winter  c xsolph - phase of year, assuming time zero is mid winter
95  c xinxx - cos ( xsolph )  c xinxx - cos ( xsolph )
96        _RL xsolph,xinxx          _RL xsolph,xinxx
97        INTEGER nForcingPeriods,Imytm,Ifprd,Ifcyc,Iftm          INTEGER nForcingPeriods,Imytm,Ifprd,Ifcyc,Iftm
98  c coefficients used to compute saturation specific humidity  c coefficients used to compute saturation specific humidity
99        DATA   ssq0,           ssq1,           ssq2          DATA   ssq0,           ssq1,           ssq2
100       &     / 3.797915 _d 0 , 7.93252 _d -6 , 2.166847 _d -3 /       &     / 3.797915 _d 0 , 7.93252 _d -6 , 2.166847 _d -3 /
101    
102        IF ( periodicExternalForcing ) THEN        IF ( periodicExternalForcing_cheap ) THEN
103    
104    
105  c the objective here is to give cheapaml a default periodic forcing  c the objective here is to give cheapaml a default periodic forcing
106  c consisting only of annually varying solar forcing, and thus Trelaxation  c consisting only of annually varying solar forcing, and thus Trelaxation
# Line 106  c keys off of solardata.  if a solar dat Line 109  c keys off of solardata.  if a solar dat
109  c assume there are files to be read and interpolated between, as is standard  c assume there are files to be read and interpolated between, as is standard
110  c for the MITGCM.  c for the MITGCM.
111    
112         IF ( SolarFile .EQ. ' '  ) THEN        IF ( SolarFile .EQ. ' '  ) THEN
113          IF (useStressOption)then        IF (useStressOption)then
114           write(*,*)' stress option is turned on.  this is not',          write(*,*)' stress option is turned on.  this is not
115       &   'consistent with the default time dependent forcing option'       & consistent with the default time dependent forcing option'
116           STOP 'ABNORMAL END: S/R CHEAPAML_FIELDS_LOAD'          stop
117          ENDIF          ENDIF
118          if ( myIter .EQ. nIter0 )then           if ( myIter .EQ. nIter0 )then
119           WRITE(*,*)           WRITE(*,*)
120       &   'S/R  Assuming Standard Annually Varying Solar Forcing'       &  'S/R  Assuming Standard Annually Varying Solar Forcing'
121          endif           endif
122          xsolph=myTime*2.d0*3.14159 _d 0/365. _d 0/86400. _d 0           xsolph=myTime*2.d0*3.14159 _d 0/365. _d 0/86400. _d 0
123          xinxx=cos(xsolph+xphaseinit+3.14159 _d 0)           xinxx=cos(xsolph+xphaseinit+3.14159 _d 0)
124          DO bj = myByLo(myThid), myByHi(myThid)           DO bj = myByLo(myThid), myByHi(myThid)
125           DO bi = myBxLo(myThid), myBxHi(myThid)            DO bi = myBxLo(myThid), myBxHi(myThid)
126            DO j=1,sNy             DO j=1,sNy
127             DO i=1,sNx              DO i=1,sNx
128              jG = myYGlobalLo-1+(bj-1)*sNy+j                    jG = myYGlobalLo-1+(bj-1)*sNy+j
129              local=225.d0+dsolms*xinxx-float((jg-1))/float((ny-1))*              local=225.d0+dsolms*xinxx-float((jg-1))/float((ny-1))*
130       &        (37.5d0-dsolmn*xinxx)       &        (37.5d0-dsolmn*xinxx)
131              Solar(i,j,bi,bj) = local                    Solar(i,j,bi,bj) = local
132                ENDDO
133             ENDDO             ENDDO
134            ENDDO            ENDDO
135           ENDDO           ENDDO
136          ENDDO           _EXCH_XY_RS(solar, mythid)
         _EXCH_XY_RL(solar, myThid)  
137  c relaxation temperature in radiative equilibrium  c relaxation temperature in radiative equilibrium
138          DO bj = myByLo(myThid), myByHi(myThid)           DO bj = myByLo(myThid), myByHi(myThid)
139           DO bi = myBxLo(myThid), myBxHi(myThid)            DO bi = myBxLo(myThid), myBxHi(myThid)
140            DO j=1,sNy             DO j=1,sNy
141             DO i=1,sNx              DO i=1,sNx
142              jG = myYGlobalLo-1+(bj-1)*sNy+j            jG = myYGlobalLo-1+(bj-1)*sNy+j
143              local=solar(i,j,bi,bj)            local=solar(i,j,bi,bj)
144              local=(2.d0*local/stefan)**(0.25d0)-Celcius2K            local=(2.d0*local/stefan)**(0.25d0)-Celsius2K
145              Tr(i,j,bi,bj) = local            TR(i,j,bi,bj) = local
146                ENDDO
147             ENDDO             ENDDO
148            ENDDO            ENDDO
149           ENDDO           ENDDO
150          ENDDO           _EXCH_XY_RS(TR, mythid)
         _EXCH_XY_RL(Tr, myThid)  
151  c default specific humidity profile to 80% relative humidity  c default specific humidity profile to 80% relative humidity
152          DO bj = myByLo(myThid), myByHi(myThid)           DO bj = myByLo(myThid), myByHi(myThid)
153           DO bi = myBxLo(myThid), myBxHi(myThid)            DO bi = myBxLo(myThid), myBxHi(myThid)
154            DO j=1,sNy             DO j=1,sNy
155             DO i=1,sNx              DO i=1,sNx
156  c           jG = myYGlobalLo-1+(bj-1)*sNy+j  c                  jG = myYGlobalLo-1+(bj-1)*sNy+j
157              local = Tr(i,j,bi,bj)+Celcius2K                    local = Tr(i,j,bi,bj)+Celsius2K
158              ssqa = ssq0*exp( lath*(ssq1-ssq2/local)) / p0                ssqa = ssq0*exp( lath*(ssq1-ssq2/local)) / p0
159              qr(i,j,bi,bj) = 0.8d0*ssqa                    qr(i,j,bi,bj) = 0.8d0*ssqa
160                ENDDO
161             ENDDO             ENDDO
162            ENDDO            ENDDO
163           ENDDO           ENDDO
164          ENDDO           _EXCH_XY_RS(qr, mythid)
         _EXCH_XY_RL(qr, myThid)  
165  c u wind field  c u wind field
166          DO bj = myByLo(myThid), myByHi(myThid)           DO bj = myByLo(myThid), myByHi(myThid)
167           DO bi = myBxLo(myThid), myBxHi(myThid)            DO bi = myBxLo(myThid), myBxHi(myThid)
168            DO j=1,sNy             DO j=1,sNy
169             DO i=1,sNx              DO i=1,sNx
170              jG = myYGlobalLo-1+(bj-1)*sNy+j                    jG = myYGlobalLo-1+(bj-1)*sNy+j
171              local=-5.d0*cos(2.d0*pi*float(jg-1)/(float(ny-1)))                    local=-5.d0*cos(2.d0*pi*float(jg-1)/(float(ny-1)))
172              uwind(i,j,bi,bj) = local                    uwind(i,j,bi,bj) = local
173                ENDDO
174             ENDDO             ENDDO
175            ENDDO            ENDDO
176           ENDDO           ENDDO
177          ENDDO           _EXCH_XY_RS(uwind, mythid)
         _EXCH_XY_RL(uwind, myThid)  
178  c v wind field  c v wind field
179          DO bj = myByLo(myThid), myByHi(myThid)           DO bj = myByLo(myThid), myByHi(myThid)
180           DO bi = myBxLo(myThid), myBxHi(myThid)            DO bi = myBxLo(myThid), myBxHi(myThid)
181            DO j=1,sNy             DO j=1,sNy
182             DO i=1,sNx              DO i=1,sNx
183              jG = myYGlobalLo-1+(bj-1)*sNy+j                    jG = myYGlobalLo-1+(bj-1)*sNy+j
184              vwind(i,j,bi,bj) = 0.d0                    vwind(i,j,bi,bj) = 0.d0
185                ENDDO
186             ENDDO             ENDDO
187            ENDDO            ENDDO
188           ENDDO           ENDDO
189          ENDDO           _EXCH_XY_RS(vwind, mythid)
190          _EXCH_XY_RL(vwind, myThid)        ELSE
   
        ELSE  
 C      else: solarFile non empty  
191    
192  C here for usual interpolative forcings  c here for usual interpolative forcings
193  C First call requires that we initialize everything to zero for safety  C First call requires that we initialize everything to zero for safety
194          IF ( myIter .EQ. nIter0 ) THEN        IF ( myIter .EQ. nIter0 ) THEN
195           DO bj = myByLo(myThid), myByHi(myThid)          DO bj = myByLo(myThid), myByHi(myThid)
196            DO bi = myBxLo(myThid), myBxHi(myThid)           DO bi = myBxLo(myThid), myBxHi(myThid)
197             DO j=1-OLy,sNy+OLy            DO j=1-OLy,sNy+OLy
198              DO i=1-OLx,sNx+OLx             DO i=1-OLx,sNx+OLx
199                trair0  (i,j,bi,bj) = 0.             trair0  (i,j,bi,bj) = 0.
200                trair1  (i,j,bi,bj) = 0.             trair1  (i,j,bi,bj) = 0.
201                qrair0  (i,j,bi,bj) = 0.             qrair0  (i,j,bi,bj) = 0.
202                qrair1  (i,j,bi,bj) = 0.             qrair1  (i,j,bi,bj) = 0.
203                solar0  (i,j,bi,bj) = 0.             solar0  (i,j,bi,bj) = 0.
204                solar1  (i,j,bi,bj) = 0.             solar1  (i,j,bi,bj) = 0.
205                uwind0  (i,j,bi,bj) = 0.             uwind0  (i,j,bi,bj) = 0.
206                uwind1  (i,j,bi,bj) = 0.             uwind1  (i,j,bi,bj) = 0.
207                vwind0  (i,j,bi,bj) = 0.             vwind0  (i,j,bi,bj) = 0.
208                vwind1  (i,j,bi,bj) = 0.             vwind1  (i,j,bi,bj) = 0.
209                ustress0(i,j,bi,bj) = 0.             ustress0(i,j,bi,bj) = 0.
210                ustress1(i,j,bi,bj) = 0.             ustress1(i,j,bi,bj) = 0.
211                vstress0(i,j,bi,bj) = 0.             vstress0(i,j,bi,bj) = 0.
212                vstress1(i,j,bi,bj) = 0.             vstress1(i,j,bi,bj) = 0.
213                wavesh0 (i,j,bi,bj) = 0.             wavesh0 (i,j,bi,bj) = 0.
214                wavesh1 (i,j,bi,bj) = 0.             wavesh1 (i,j,bi,bj) = 0.
215                wavesp0 (i,j,bi,bj) = 0.             wavesp0 (i,j,bi,bj) = 0.
216                wavesp1 (i,j,bi,bj) = 0.             wavesp1 (i,j,bi,bj) = 0.
217                rair0   (i,j,bi,bj) = 0.             rair0   (i,j,bi,bj) = 0.
218                rair1   (i,j,bi,bj) = 0.             rair1   (i,j,bi,bj) = 0.
             ENDDO  
219             ENDDO             ENDDO
220            ENDDO            ENDDO
221           ENDDO           ENDDO
222          ENDIF          ENDDO      
223          ENDIF
224    
225  C Now calculate whether it is time to update the forcing arrays  C Now calculate whether it is time to update the forcing arrays
226          rdt=1. _d 0 / deltaTclock        rdt=1. _d 0 / deltaTclock
227          nForcingPeriods=        nForcingPeriods=
228       &   int(externForcingCycle/externForcingPeriod+0.5)       &  int(externForcingCycle_cheap/externForcingPeriod_cheap+0.5)
229          Imytm=int(myTime*rdt+0.5)        Imytm=int(myTime*rdt+0.5)
230          Ifprd=int(externForcingPeriod*rdt+0.5)        Ifprd=int(externForcingPeriod_cheap*rdt+0.5)
231          Ifcyc=int(externForcingCycle*rdt+0.5)        Ifcyc=int(externForcingCycle_cheap*rdt+0.5)
232          Iftm=mod( Imytm+Ifcyc-Ifprd/2 ,Ifcyc)        Iftm=mod( Imytm+Ifcyc-Ifprd/2 ,Ifcyc)
233    
234          intime0=int(Iftm/Ifprd)        intime0=int(Iftm/Ifprd)
235          intime1=mod(intime0+1,nForcingPeriods)        intime1=mod(intime0+1,nForcingPeriods)
236  c       aWght=float( Iftm-Ifprd*intime0 )/float( Ifprd )  c     aWght=float( Iftm-Ifprd*intime0 )/float( Ifprd )
237          aWght=dfloat( Iftm-Ifprd*intime0 )/dfloat( Ifprd )        aWght=dfloat( Iftm-Ifprd*intime0 )/dfloat( Ifprd )
238          bWght=1.-aWght        bWght=1.-aWght
239    
240          intime0=intime0+1        intime0=intime0+1
241          intime1=intime1+1        intime1=intime1+1
242    
243          IF (        IF (
244       &       Iftm-Ifprd*(intime0-1) .EQ. 0       &  Iftm-Ifprd*(intime0-1) .EQ. 0
245       &      .OR. myIter .EQ. nIter0       &  .OR. myIter .EQ. nIter0
246       &     ) THEN       & ) THEN
247    
248  C      If the above condition is met then we need to read in  C      If the above condition is met then we need to read in
249  C      data for the period ahead and the period behind myTime.  C      data for the period ahead and the period behind myTime.
250           WRITE(*,*)         WRITE(*,*)
251       &     'S/R CHEAPAML_FIELDS_LOAD'       &  'S/R CHEAPAML_FIELDS_LOAD'
252           IF ( SolarFile .NE. ' '  ) THEN        IF ( SolarFile .NE. ' '  ) THEN
253            CALL READ_REC_XY_RL( SolarFile,solar0,intime0,         CALL READ_REC_XY_RL( SolarFile,solar0,intime0,
254       &         myIter,myThid )       &       myIter,myThid )
255            CALL READ_REC_XY_RL( SolarFile,solar1,intime1,         CALL READ_REC_XY_RL( SolarFile,solar1,intime1,
256       &         myIter,myThid )       &       myIter,myThid )
257           ENDIF        ENDIF
258           IF ( TrFile .NE. ' '  ) THEN        IF ( TrFile .NE. ' '  ) THEN
259            CALL READ_REC_XY_RL( TrFile,trair0,intime0,         CALL READ_REC_XY_RL( TRFile,trair0,intime0,
260       &         myIter,myThid )       &       myIter,myThid )
261            CALL READ_REC_XY_RL( TrFile,trair1,intime1,         CALL READ_REC_XY_RL( TRFile,trair1,intime1,
262       &         myIter,myThid )       &       myIter,myThid )
263           ENDIF        ENDIF
264           IF ( QrFile .NE. ' '  ) THEN        IF ( QrFile .NE. ' '  ) THEN
265            CALL READ_REC_XY_RL( QrFile,qrair0,intime0,         CALL READ_REC_XY_RL( QrFile,qrair0,intime0,
266       &         myIter,myThid )       &       myIter,myThid )
267            CALL READ_REC_XY_RL( QrFile,qrair1,intime1,         CALL READ_REC_XY_RL( QrFile,qrair1,intime1,
268       &         myIter,myThid )       &       myIter,myThid )
269           ENDIF        ENDIF
270           IF ( UWindFile .NE. ' '  ) THEN        IF ( UWindFile .NE. ' '  ) THEN
271            CALL READ_REC_XY_RL( UWindFile,uwind0,intime0,         CALL READ_REC_XY_RL( UWindFile,uwind0,intime0,
272       &         myIter,myThid )       &       myIter,myThid )
273            CALL READ_REC_XY_RL( UWindFile,uwind1,intime1,         CALL READ_REC_XY_RL( UWindFile,uwind1,intime1,
274       &         myIter,myThid )       &       myIter,myThid )
275           ENDIF        ENDIF
276           IF ( VWindFile .NE. ' '  ) THEN        IF ( VWindFile .NE. ' '  ) THEN
277            CALL READ_REC_XY_RL( VWindFile,vwind0,intime0,         CALL READ_REC_XY_RL( VWindFile,vwind0,intime0,
278       &         myIter,myThid )       &       myIter,myThid )
279            CALL READ_REC_XY_RL( VWindFile,vwind1,intime1,         CALL READ_REC_XY_RL( VWindFile,vwind1,intime1,
280       &         myIter,myThid )       &       myIter,myThid )
281           ENDIF        ENDIF
282           IF(useStressOption)THEN        IF(useStressOption)THEN
283            IF ( UStressFile .NE. ' '  ) THEN        IF ( UStressFile .NE. ' '  ) THEN
284             CALL READ_REC_XY_RL( UStressFile,ustress0,intime0,         CALL READ_REC_XY_RL( UStressFile,ustress0,intime0,
285       &          myIter,myThid )       &       myIter,myThid )
286             CALL READ_REC_XY_RL( UStressFile,ustress1,intime1,         CALL READ_REC_XY_RL( UStressFile,ustress1,intime1,
287       &          myIter,myThid )       &       myIter,myThid )
288            ENDIF        ENDIF
289            IF ( VStressFile .NE. ' '  ) THEN        IF ( VStressFile .NE. ' '  ) THEN
290             CALL READ_REC_XY_RL( VStressFile,vstress0,intime0,         CALL READ_REC_XY_RL( VStressFile,vstress0,intime0,
291       &          myIter,myThid )       &       myIter,myThid )
292             CALL READ_REC_XY_RL( VStressFile,vstress1,intime1,         CALL READ_REC_XY_RL( VStressFile,vstress1,intime1,
293       &          myIter,myThid )       &       myIter,myThid )
294            ENDIF          ENDIF
295           ENDIF        ENDIF
296           IF(useRelativeHumidity)THEN        IF ( FluxFormula.eq.'COARE3') THEN
297            IF ( QrrelFile .NE. ' '  ) THEN        IF ( WaveHFile .NE. ' '  ) THEN
298             CALL READ_REC_XY_RL( QrrelFile,rair0,intime0,         CALL READ_REC_XY_RL( WaveHFile,wavesh0,intime0,
299       &          myIter,myThid )       &       myIter,myThid )
300             CALL READ_REC_XY_RL( QrrelFile,rair1,intime1,         CALL READ_REC_XY_RL( WaveHFile,wavesh1,intime1,
301       &          myIter,myThid )       &       myIter,myThid )
302            ENDIF        ENDIF
303  C     This subroutine is in cheapaml_init_varia        IF ( WavePFile .NE. ' '  ) THEN
304            CALL CHEAPAML_CONVERT_HUM(rair0, trair0, qrair0,myThid )         CALL READ_REC_XY_RL( WavePFile,wavesp0,intime0,
305            CALL CHEAPAML_CONVERT_HUM(rair1, trair1, qrair1,myThid )       &       myIter,myThid )
306           ENDIF         CALL READ_REC_XY_RL( WavePFile,wavesp1,intime1,
307           IF ( FluxFormula.EQ.'COARE3') THEN       &       myIter,myThid )
308            IF ( WaveHFile .NE. ' '  ) THEN        ENDIF
309             CALL READ_REC_XY_RL( WaveHFile,wavesh0,intime0,          ENDIF
      &          myIter,myThid )  
            CALL READ_REC_XY_RL( WaveHFile,wavesh1,intime1,  
      &          myIter,myThid )  
           ENDIF  
           IF ( WavePFile .NE. ' '  ) THEN  
            CALL READ_REC_XY_RL( WavePFile,wavesp0,intime0,  
      &          myIter,myThid )  
            CALL READ_REC_XY_RL( WavePFile,wavesp1,intime1,  
      &          myIter,myThid )  
           ENDIF  
          ENDIF  
   
          _EXCH_XY_RL( trair0 , myThid )  
          _EXCH_XY_RL( qrair0 , myThid )  
          _EXCH_XY_RL( solar0 , myThid )  
          _EXCH_XY_RL( uwind0 , myThid )  
          _EXCH_XY_RL( vwind0 , myThid )  
          _EXCH_XY_RL( trair1 , myThid )  
          _EXCH_XY_RL( qrair1 , myThid )  
          _EXCH_XY_RL( solar1 , myThid )  
          _EXCH_XY_RL( uwind1 , myThid )  
          _EXCH_XY_RL( vwind1 , myThid )  
          IF(useStressOption)THEN  
           _EXCH_XY_RL( ustress0 , myThid )  
           _EXCH_XY_RL( vstress0 , myThid )  
           _EXCH_XY_RL( ustress1 , myThid )  
           _EXCH_XY_RL( vstress1 , myThid )  
          ENDIF  
          IF(FluxFormula.EQ.'COARE3') THEN  
           _EXCH_XY_RL( wavesp0 , myThid )  
           _EXCH_XY_RL( wavesp1 , myThid )  
           _EXCH_XY_RL( wavesh0 , myThid )  
           _EXCH_XY_RL( wavesh1 , myThid )  
          ENDIF  
310    
311  C       end of loading new fields block         _EXCH_XY_RL(trair0 , myThid )
312           _EXCH_XY_RL(qrair0 , myThid )
313           _EXCH_XY_RL(solar0 , myThid )
314           _EXCH_XY_RL(uwind0 , myThid )
315           _EXCH_XY_RL(vwind0 , myThid )
316           _EXCH_XY_RL(trair1 , myThid )
317           _EXCH_XY_RL(qrair1 , myThid )
318           _EXCH_XY_RL(solar1 , myThid )
319           _EXCH_XY_RL(uwind1 , myThid )
320           _EXCH_XY_RL(vwind1 , myThid )
321           IF(useStressOption)THEN
322           _EXCH_XY_RL(uwind0 , myThid )
323           _EXCH_XY_RL(vwind0 , myThid )
324           _EXCH_XY_RL(uwind1 , myThid )
325           _EXCH_XY_RL(vwind1 , myThid )
326            ENDIF
327            IF(FluxFormula.EQ.'COARE3') THEN
328           _EXCH_XY_RS(wavesp0 , myThid )
329           _EXCH_XY_RS(wavesp1 , myThid )
330           _EXCH_XY_RS(wavesh0 , myThid )
331           _EXCH_XY_RS(wavesh1 , myThid )
332          ENDIF          ENDIF
333    
334  C--   Interpolate Tr, Qr, Solar  C     end of loading new fields block
335          DO bj = myByLo(myThid), myByHi(myThid)        ENDIF
336           DO bi = myBxLo(myThid), myBxHi(myThid)  
337            DO j=1,sNy  C--   Interpolate TR, QR, SOLAR
338             DO i=1,sNx        DO bj = myByLo(myThid), myByHi(myThid)
339              Tr(i,j,bi,bj)   = bWght*trair0(i,j,bi,bj)         DO bi = myBxLo(myThid), myBxHi(myThid)
340       &                       +aWght*trair1(i,j,bi,bj)   !+273.15          DO j=1-OLy,sNy+OLy
341              qr(i,j,bi,bj)   = bWght*qrair0(i,j,bi,bj)           DO i=1-OLx,sNx+OLx
342       &                       +aWght*qrair1(i,j,bi,bj)            TR(i,j,bi,bj)      = bWght*trair0(i,j,bi,bj)
343              uwind(i,j,bi,bj)= bWght*uwind0(i,j,bi,bj)       &                        +aWght*trair1(i,j,bi,bj)   !+273.15
344       &                       +aWght*uwind1(i,j,bi,bj)            qr(i,j,bi,bj)      = bWght*qrair0(i,j,bi,bj)
345              vwind(i,j,bi,bj)= bWght*vwind0(i,j,bi,bj)       &                        +aWght*qrair1(i,j,bi,bj)
346       &                       +aWght*vwind1(i,j,bi,bj)            uwind(i,j,bi,bj)   = bWght*uwind0(i,j,bi,bj)
347              solar(i,j,bi,bj)= bWght*solar0(i,j,bi,bj)       &                        +aWght*uwind1(i,j,bi,bj)
348       &                       +aWght*solar1(i,j,bi,bj)            vwind(i,j,bi,bj)   = bWght*vwind0(i,j,bi,bj)
349              IF(useStressOption)THEN       &                        +aWght*vwind1(i,j,bi,bj)
350               ustress(i,j,bi,bj)= bWght*ustress0(i,j,bi,bj)            solar(i,j,bi,bj)   = bWght*solar0(i,j,bi,bj)
351         &                        +aWght*solar1(i,j,bi,bj)
352           IF(useStressOption)THEN
353              ustress(i,j,bi,bj)   = bWght*ustress0(i,j,bi,bj)
354       &                          +aWght*ustress1(i,j,bi,bj)       &                          +aWght*ustress1(i,j,bi,bj)
355               vstress(i,j,bi,bj)= bWght*vstress0(i,j,bi,bj)            vstress(i,j,bi,bj)   = bWght*vstress0(i,j,bi,bj)
356       &                          +aWght*vstress1(i,j,bi,bj)       &                          +aWght*vstress1(i,j,bi,bj)
357              ENDIF         ENDIF
358              IF(FluxFormula.EQ.'COARE3')THEN          IF(FluxFormula.eq.'COARE3')THEN
359               IF(WaveHFile.NE.' ')THEN          IF(WaveHFile.ne.' ')THEN
360                wavesh(i,j,bi,bj) = bWght*wavesh0(i,j,bi,bj)            wavesh(i,j,bi,bj)   = bWght*wavesh0(i,j,bi,bj)
361       &                           +aWght*wavesh1(i,j,bi,bj)       &                         +aWght*wavesh1(i,j,bi,bj)
362               ENDIF          ENDIF
363               IF(WavePFile.NE.' ')THEN          IF(WavePFile.ne.' ')THEN
364                wavesp(i,j,bi,bj) = bWght*wavesp0(i,j,bi,bj)            wavesp(i,j,bi,bj)   = bWght*wavesp0(i,j,bi,bj)
365       &                           +aWght*wavesp1(i,j,bi,bj)       &                         +aWght*wavesp1(i,j,bi,bj)
366               ENDIF          ENDIF
367              ELSE          ELSE
368               u=uwind(i,j,bi,bj)**2+vwind(i,j,bi,bj)**2            u=uwind(i,j,bi,bj)**2+vwind(i,j,bi,bj)**2
369               u=dsqrt(u)            u=dsqrt(u)
370               wavesp(i,j,bi,bj)=0.729 _d 0 * u            wavesp(i,j,bi,bj)=0.729 _d 0 * u  
371               wavesh(i,j,bi,bj)=0.018 _d 0 * u*u*(1. _d 0 + .015 _d 0 *u)            wavesh(i,j,bi,bj)=0.018 _d 0 * u*u*(1. _d 0 + .015 _d 0 *u)
372              ENDIF            ENDIF
            ENDDO  
           ENDDO  
373           ENDDO           ENDDO
374          ENDDO          ENDDO
375           ENDDO
376          ENDDO
377    
378  C end if solarFile is empty  c end if solarFile is empty      
379         ENDIF        ENDIF
380    
381    c end of periodic forcing options, on to steady option
382    
 C end of periodic forcing options, on to steady option  
383        ELSE        ELSE
384    
385         IF ( myIter .EQ. nIter0 ) THEN         IF ( myIter .EQ. nIter0 ) THEN
386          IF ( SolarFile .NE. ' '  ) THEN          IF ( SolarFile .NE. ' '  ) THEN
387           CALL READ_FLD_XY_RL( SolarFile,' ',solar,0,myThid )           CALL READ_FLD_XY_RS( SolarFile,' ',solar,0,myThid )
388          ELSE          ELSE
389           DO bj = myByLo(myThid), myByHi(myThid)           DO bj = myByLo(myThid), myByHi(myThid)
390            DO bi = myBxLo(myThid), myBxHi(myThid)            DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 407  C end of periodic forcing options, on to Line 398  C end of periodic forcing options, on to
398            ENDDO            ENDDO
399           ENDDO           ENDDO
400          ENDIF          ENDIF
401          _EXCH_XY_RL(solar, myThid)          _EXCH_XY_RS(solar, mythid)
402          IF ( TrFile .NE. ' '  ) THEN          IF ( TrFile .NE. ' '  ) THEN
403           CALL READ_FLD_XY_RL( TrFile,' ',tr,0,myThid )           CALL READ_FLD_XY_RL( TrFile,' ',tr,0,myThid )
404          ELSE          ELSE
# Line 418  C end of periodic forcing options, on to Line 409  C end of periodic forcing options, on to
409              jG = myYGlobalLo-1+(bj-1)*sNy+j              jG = myYGlobalLo-1+(bj-1)*sNy+j
410              local=solar(i,j,bi,bj)              local=solar(i,j,bi,bj)
411              local=(2.d0*local/stefan)**(0.25d0)-273.16              local=(2.d0*local/stefan)**(0.25d0)-273.16
412              Tr(i,j,bi,bj) = local              TR(i,j,bi,bj) = local
413              ENDDO              ENDDO
414             ENDDO             ENDDO
415            ENDDO            ENDDO
416           ENDDO           ENDDO
417          ENDIF          ENDIF
418          _EXCH_XY_RL(Tr, myThid)          
419            _EXCH_XY_RL(TR, mythid)
420  C do specific humidity  c do specific humidity
421          IF ( QrFile .NE. ' '.AND..NOT. useRelativeHumidity  ) THEN          IF ( QrFile .NE. ' ') THEN
422           CALL READ_FLD_XY_RL( QrFile,' ',qr,0,myThid )           CALL READ_FLD_XY_RL( QrFile,' ',qr,0,myThid )      
         ELSEIF ( QrrelFile .NE. ' '.AND.useRelativeHumidity) THEN  
          CALL READ_FLD_XY_RL( QrrelFile,' ',rref,0,myThid )  
          CALL CHEAPAML_CONVERT_HUM(rref, Tr, qr,myThid )  
   
423          ELSE          ELSE
424  C default specific humidity profile to 80% relative humidity  c default specific humidity profile to 80% relative humidity
425           DO bj = myByLo(myThid), myByHi(myThid)           DO bj = myByLo(myThid), myByHi(myThid)
426            DO bi = myBxLo(myThid), myBxHi(myThid)            DO bi = myBxLo(myThid), myBxHi(myThid)
427             DO j=1,sNy             DO j=1,sNy
# Line 448  c                  jG = myYGlobalLo-1+(b Line 435  c                  jG = myYGlobalLo-1+(b
435            ENDDO            ENDDO
436           ENDDO           ENDDO
437          ENDIF          ENDIF
438          _EXCH_XY_RL(qr, myThid)          _EXCH_XY_RL(qr, mythid)
439          IF ( UWindFile .NE. ' '  ) THEN          IF ( UWindFile .NE. ' '  ) THEN
440           CALL READ_FLD_XY_RL( UWindFile,' ',uwind,0,myThid )           CALL READ_FLD_XY_RL( UWindFile,' ',uwind,0,myThid )
441          ELSE          ELSE
# Line 456  c                  jG = myYGlobalLo-1+(b Line 443  c                  jG = myYGlobalLo-1+(b
443            DO bi = myBxLo(myThid), myBxHi(myThid)            DO bi = myBxLo(myThid), myBxHi(myThid)
444             DO j=1,sNy             DO j=1,sNy
445              DO i=1,sNx              DO i=1,sNx
446                jG = myYGlobalLo-1+(bj-1)*sNy+j                    jG = myYGlobalLo-1+(bj-1)*sNy+j
447  c mod for debug  c mod for debug
448  c to return to original code, uncomment following line  c to return to original code, uncomment following line
449  c comment out 2nd line  c comment out 2nd line
450                local=-5.d0*cos(2.d0*pi*float(jg-1)/(float(ny-1)))                    local=-5.d0*cos(2.d0*pi*float(jg-1)/(float(ny-1)))
451  c             local=0.d0*cos(2.d0*pi*float(jg-1)/(float(ny-1)))  c                 local=0.d0*cos(2.d0*pi*float(jg-1)/(float(ny-1)))
452                uwind(i,j,bi,bj) = local                    uwind(i,j,bi,bj) = local
453              ENDDO              ENDDO
454             ENDDO             ENDDO
455            ENDDO            ENDDO
456           ENDDO           ENDDO
457          ENDIF          ENDIF
458           _EXCH_XY_RL(uwind, myThid)           _EXCH_XY_RL(uwind, mythid)
459          IF ( VWindFile .NE. ' '  ) THEN          IF ( VWindFile .NE. ' '  ) THEN
460           CALL READ_FLD_XY_RL( VWindFile,' ',vwind,0,myThid )           CALL READ_FLD_XY_RL( VWindFile,' ',vwind,0,myThid )
461          ELSE          ELSE
# Line 476  c             local=0.d0*cos(2.d0*pi*flo Line 463  c             local=0.d0*cos(2.d0*pi*flo
463            DO bi = myBxLo(myThid), myBxHi(myThid)            DO bi = myBxLo(myThid), myBxHi(myThid)
464             DO j=1,sNy             DO j=1,sNy
465              DO i=1,sNx              DO i=1,sNx
466                jG = myYGlobalLo-1+(bj-1)*sNy+j                    jG = myYGlobalLo-1+(bj-1)*sNy+j
467                vwind(i,j,bi,bj) = 0.d0                    vwind(i,j,bi,bj) = 0.d0
468              ENDDO              ENDDO
469             ENDDO             ENDDO
470            ENDDO            ENDDO
471           ENDDO           ENDDO
472          ENDIF          ENDIF
473          _EXCH_XY_RL(vwind, myThid)           _EXCH_XY_RL(vwind, mythid)
474          IF(useStressOption)THEN         IF(useStressOption)THEN
475           IF ( UStressFile .NE. ' '  ) THEN          IF ( UStressFile .NE. ' '  ) THEN
476            CALL READ_FLD_XY_RL( UStressFile,' ',ustress,0,myThid )           CALL READ_FLD_XY_RL( UStressFile,' ',ustress,0,myThid )
477           ELSE          ELSE
478            write(*,*)' U Stress File absent with stress option'          write(*,*)' U Stress File absent with stress option'
479            STOP 'ABNORMAL END: S/R CHEAPAML_FIELDS_LOAD'          stop
480           ENDIF          ENDIF
481           IF ( VStressFile .NE. ' '  ) THEN          IF ( VStressFile .NE. ' '  ) THEN
482            CALL READ_FLD_XY_RL( VStressFile,' ',vstress,0,myThid )           CALL READ_FLD_XY_RS( VStressFile,' ',vstress,0,myThid )
483           ELSE          ELSE
484            write(*,*)' V Stress File absent with stress option'          write(*,*)' V Stress File absent with stress option'
485            STOP 'ABNORMAL END: S/R CHEAPAML_FIELDS_LOAD'          stop
486           ENDIF          ENDIF
487           _EXCH_XY_RL(ustress, myThid)           _EXCH_XY_RL(ustress, mythid)
488           _EXCH_XY_RL(vstress, myThid)           _EXCH_XY_RL(vstress, mythid)
489          ENDIF          ENDIF
490          IF (FluxFormula.EQ.'COARE3')THEN          IF (FluxFormula.eq.'COARE3')THEN
491           IF (WaveHFile.NE.' ')THEN                  IF (WaveHFile.NE.' ')THEN
492            CALL READ_FLD_XY_RL( WaveHFile,' ',wavesh,0,myThid )           CALL READ_FLD_XY_RL( WaveHFile,' ',wavesh,0,myThid )
493           ENDIF          ENDIF
494           IF (WavePFile.NE.' ')THEN                  IF (WavePFile.NE.' ')THEN
495            CALL READ_FLD_XY_RL( WavePFile,' ',wavesp,0,myThid )           CALL READ_FLD_XY_RL( WavePFile,' ',wavesp,0,myThid )
496           ELSE          ELSE
497            DO bj = myByLo(myThid), myByHi(myThid)           DO bj = myByLo(myThid), myByHi(myThid)
498             DO bi = myBxLo(myThid), myBxHi(myThid)            DO bi = myBxLo(myThid), myBxHi(myThid)
499              DO j=1,sNy             DO j=1,sNy
500               DO i=1,sNx              DO i=1,sNx
501                u=uwind(i,j,bi,bj)**2+vwind(i,j,bi,bj)**2                  u=uwind(i,j,bi,bj)**2+vwind(i,j,bi,bj)**2
502                u=dsqrt(u)                  u=dsqrt(u)
503                wavesp(i,j,bi,bj)=0.729 _d 0*u            wavesp(i,j,bi,bj)=0.729 _d 0 * u  
504                wavesh(i,j,bi,bj)=0.018 _d 0*u*u*(1. _d 0 + .015 _d 0 *u)            wavesh(i,j,bi,bj)=0.018 _d 0 * u*u*(1. _d 0 + .015 _d 0 *u)
              ENDDO  
505              ENDDO              ENDDO
506             ENDDO             ENDDO
507            ENDDO            ENDDO
508           ENDIF           ENDDO
509           _EXCH_XY_RL(wavesp, myThid)         ENDIF
510           _EXCH_XY_RL(wavesh, myThid)          _EXCH_XY_RL(wavesp, mythid)
511          ENDIF          _EXCH_XY_RL(wavesh, mythid)
512          ENDIF
 C end if myIter = nIter0  
513         ENDIF         ENDIF
514    
515    
516  C endif for Steady Option  C endif for Steady Option
517        ENDIF        ENDIF
518    
519  C fill in outer edges  C fill in outer edges
520        DO bj = myByLo(myThid), myByHi(myThid)  
521          DO bi = myBxLo(myThid), myBxHi(myThid)           DO bj = myByLo(myThid), myByHi(myThid)
522            DO j=1-oly,sny+oly            DO bi = myBxLo(myThid), myBxHi(myThid)
523              jG = myYGlobalLo-1+(bj-1)*sNy+j                  do j=1-oly,sny+oly
524              DO i=1-olx,snx+olx                    jG = myYGlobalLo-1+(bj-1)*sNy+j
525                iG=myXGlobalLo-1+(bi-1)*sNx+i                  do i=1-olx,snx+olx
526                if(iG.lt.1)then                    iG=myXGlobalLo-1+(bi-1)*sNx+i
527                  Tr(i,j,bi,bj)=Tr(1,j,bi,bj)              if(iG.lt.1)then
528                  qr(i,j,bi,bj)=qr(1,j,bi,bj)              Tr(i,j,bi,bj)=Tr(1,j,bi,bj)
529                  uwind(i,j,bi,bj)=uwind(1,j,bi,bj)              qr(i,j,bi,bj)=qr(1,j,bi,bj)
530                  vwind(i,j,bi,bj)=vwind(1,j,bi,bj)              uwind(i,j,bi,bj)=uwind(1,j,bi,bj)
531                  Solar(i,j,bi,bj)=Solar(1,j,bi,bj)              vwind(i,j,bi,bj)=vwind(1,j,bi,bj)
532                Solar(i,j,bi,bj)=Solar(1,j,bi,bj)
533                
534                  if(UseStressOption)then                  if(UseStressOption)then
535                    ustress(i,j,bi,bj)=ustress(1,j,bi,bj)                  ustress(i,j,bi,bj)=ustress(1,j,bi,bj)
536                    vstress(i,j,bi,bj)=vstress(1,j,bi,bj)                  vstress(i,j,bi,bj)=vstress(1,j,bi,bj)
537                  endif                  endif
538                  if(FluxFormula.EQ.'COARE3')then  
539                    wavesp(i,j,bi,bj)=wavesp(1,j,bi,bj)                  if(FluxFormula.eq.'COARE3')then
540                    wavesh(i,j,bi,bj)=wavesh(1,j,bi,bj)                  wavesp(i,j,bi,bj)=wavesp(1,j,bi,bj)
541                    wavesh(i,j,bi,bj)=wavesh(1,j,bi,bj)
542                  endif                  endif
543                elseif(iG.gt.Nx)then  
544                  Tr(i,j,bi,bj)=Tr(sNx,j,bi,bj)              elseif(iG.gt.Nx)then
545                  qr(i,j,bi,bj)=qr(sNx,j,bi,bj)              Tr(i,j,bi,bj)=Tr(sNx,j,bi,bj)
546                  uwind(i,j,bi,bj)=uwind(sNx,j,bi,bj)              qr(i,j,bi,bj)=qr(sNx,j,bi,bj)
547                  vwind(i,j,bi,bj)=vwind(sNx,j,bi,bj)              uwind(i,j,bi,bj)=uwind(sNx,j,bi,bj)
548                  Solar(i,j,bi,bj)=Solar(sNx,j,bi,bj)              vwind(i,j,bi,bj)=vwind(sNx,j,bi,bj)
549                Solar(i,j,bi,bj)=Solar(sNx,j,bi,bj)
550    
551                  if(UseStressOption)then                  if(UseStressOption)then
552                    ustress(i,j,bi,bj)=ustress(sNx,j,bi,bj)                  ustress(i,j,bi,bj)=ustress(sNx,j,bi,bj)
553                    vstress(i,j,bi,bj)=vstress(sNx,j,bi,bj)                  vstress(i,j,bi,bj)=vstress(sNx,j,bi,bj)
554                  endif                  endif
555                  if(FluxFormula.EQ.'COARE3')then  
556                    wavesp(i,j,bi,bj)=wavesp(sNx,j,bi,bj)                  if(FluxFormula.eq.'COARE3')then
557                    wavesh(i,j,bi,bj)=wavesh(sNx,j,bi,bj)                  wavesp(i,j,bi,bj)=wavesp(sNx,j,bi,bj)
558                    wavesh(i,j,bi,bj)=wavesh(sNx,j,bi,bj)
559                  endif                  endif
560                elseif(jG.lt.1)then  
561                  Tr(i,j,bi,bj)=Tr(i,1,bi,bj)              elseif(jG.lt.1)then
562                  qr(i,j,bi,bj)=qr(i,1,bi,bj)              Tr(i,j,bi,bj)=Tr(i,1,bi,bj)
563                  uwind(i,j,bi,bj)=uwind(i,1,bi,bj)              qr(i,j,bi,bj)=qr(i,1,bi,bj)
564                  vwind(i,j,bi,bj)=vwind(i,1,bi,bj)              uwind(i,j,bi,bj)=uwind(i,1,bi,bj)
565                  Solar(i,j,bi,bj)=Solar(i,1,bi,bj)              vwind(i,j,bi,bj)=vwind(i,1,bi,bj)
566                Solar(i,j,bi,bj)=Solar(i,1,bi,bj)
567    
568                  if(UseStressOption)then                  if(UseStressOption)then
569                    ustress(i,j,bi,bj)=ustress(i,1,bi,bj)                  ustress(i,j,bi,bj)=ustress(i,1,bi,bj)
570                    vstress(i,j,bi,bj)=vstress(i,1,bi,bj)                  vstress(i,j,bi,bj)=vstress(i,1,bi,bj)
571                  endif                  endif
572                  if(FluxFormula.EQ.'COARE3')then  
573                    wavesp(i,j,bi,bj)=wavesp(i,1,bi,bj)                  if(FluxFormula.eq.'COARE3')then
574                    wavesh(i,j,bi,bj)=wavesh(i,1,bi,bj)                  wavesp(i,j,bi,bj)=wavesp(i,1,bi,bj)
575                    wavesh(i,j,bi,bj)=wavesh(i,1,bi,bj)
576                  endif                  endif
577                elseif(jG.gt.sNy)then  
578                  Tr(i,j,bi,bj)=Tr(i,sNy,bi,bj)              elseif(jG.gt.Ny)then
579                  qr(i,j,bi,bj)=qr(i,sNy,bi,bj)              Tr(i,j,bi,bj)=Tr(i,sNy,bi,bj)
580                  uwind(i,j,bi,bj)=uwind(i,sNy,bi,bj)              qr(i,j,bi,bj)=qr(i,sNy,bi,bj)
581                  vwind(i,j,bi,bj)=vwind(i,sNy,bi,bj)              uwind(i,j,bi,bj)=uwind(i,sNy,bi,bj)
582                  Solar(i,j,bi,bj)=Solar(i,sNy,bi,bj)              vwind(i,j,bi,bj)=vwind(i,sNy,bi,bj)
583                Solar(i,j,bi,bj)=Solar(i,sNy,bi,bj)
584    
585                  if(UseStressOption)then                  if(UseStressOption)then
586                    ustress(i,j,bi,bj)=ustress(i,sNy,bi,bj)                  ustress(i,j,bi,bj)=ustress(i,sNy,bi,bj)
587                    vstress(i,j,bi,bj)=vstress(i,sNy,bi,bj)                  vstress(i,j,bi,bj)=vstress(i,sNy,bi,bj)
588                  endif                  endif
589                  if(FluxFormula.EQ.'COARE3')then  
590                    wavesp(i,j,bi,bj)=wavesp(i,sNy,bi,bj)                  if(FluxFormula.eq.'COARE3')then
591                    wavesh(i,j,bi,bj)=wavesh(i,sNy,bi,bj)                  wavesp(i,j,bi,bj)=wavesp(i,sNy,bi,bj)
592                    wavesh(i,j,bi,bj)=wavesh(i,sNy,bi,bj)
593                  endif                  endif
594                endif  
595                endif
596              ENDDO              ENDDO
597               ENDDO
598            ENDDO            ENDDO
599          ENDDO           ENDDO
600        ENDDO         RETURN
   
       RETURN  
601        END        END
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22