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