/[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.2 by jmc, Tue Sep 9 18:38:46 2008 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         &         vwind0,   vwind1,
61        _RS  trair0    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)       &         ustress0, ustress1,
62        _RS  trair1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)       &         vstress0, vstress1,
63        _RS  qrair0    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)       &         wavesh0,  wavesh1,
64        _RS  qrair1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)       &         wavesp0,  wavesp1,
65        _RS  Solar0    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)       &         rair0,    rair1
66        _RS  Solar1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  
67        _RS  uwind0    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  trair0    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
68        _RS  uwind1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  trair1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
69        _RS  vwind0    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  qrair0    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
70        _RS  vwind1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  qrair1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
71          _RL  Solar0    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
72          _RL  Solar1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
73          _RL  uwind0    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
74          _RL  uwind1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
75          _RL  vwind0    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
76          _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,q        _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 80  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 96  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=1,nSy          DO bj = myByLo(myThid), myByHi(myThid)
122           DO bi=1,nSx           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
129                    Solar(i,j,bi,bj) = local             ENDDO
130                 ENDDO            ENDDO
             ENDDO  
131           ENDDO           ENDDO
132        ENDDO          ENDDO
133         _EXCH_XY_RS(solar, mythid)          _EXCH_XY_RL(solar, myThid)
134  c relaxation temperature in radiative equilibrium  c relaxation temperature in radiative equilibrium
135        DO bj=1,nSy          DO bj = myByLo(myThid), myByHi(myThid)
136           DO bi=1,nSx           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              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
143            local=local+bump             ENDDO
144            TR(i,j,bi,bj) = local            ENDDO
                ENDDO  
             ENDDO  
145           ENDDO           ENDDO
146        ENDDO          ENDDO
147        _EXCH_XY_RS(TR, mythid)          _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=1,nSy          DO bj = myByLo(myThid), myByHi(myThid)
150           DO bi=1,nSx           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
157                 ENDDO             ENDDO
158              ENDDO            ENDDO
159           ENDDO           ENDDO
160        ENDDO          ENDDO
161        _EXCH_XY_RS(qr, mythid)          _EXCH_XY_RL(qr, myThid)
162  c u wind field  c u wind field
163        DO bj=1,nSy          DO bj = myByLo(myThid), myByHi(myThid)
164           DO bi=1,nSx           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
170                 ENDDO             ENDDO
171              ENDDO            ENDDO
172           ENDDO           ENDDO
173        ENDDO          ENDDO
174        _EXCH_XY_RS(uwind, mythid)          _EXCH_XY_RL(uwind, myThid)
175  c v wind field  c v wind field
176        DO bj=1,nSy          DO bj = myByLo(myThid), myByHi(myThid)
177           DO bi=1,nSx           DO bi = myBxLo(myThid), myBxHi(myThid)
178              DO j=1,sNy            DO j=1,sNy
179                 DO i=1,sNx             DO i=1,sNx
180                    jG = myYGlobalLo-1+(bj-1)*sNy+j              jG = myYGlobalLo-1+(bj-1)*sNy+j
181                    vwind(i,j,bi,bj) = 0.d0              vwind(i,j,bi,bj) = 0.d0
182                 ENDDO             ENDDO
183              ENDDO            ENDDO
184           ENDDO           ENDDO
185        ENDDO          ENDDO
186        _EXCH_XY_RS(vwind, mythid)          _EXCH_XY_RL(vwind, myThid)
187        ELSE  
188           ELSE
189    C      else: solarFile non empty
190    
191  c here for usual interpolative forcings  C here for usual interpolative forcings
192  C First call requires that we initialize everything to zero for safety  C First call requires that we initialize everything to zero for safety
193        IF ( myIter .EQ. nIter0 ) THEN          IF ( myIter .EQ. nIter0 ) THEN
194         CALL LEF_ZERO( trair0 ,myThid )           DO bj = myByLo(myThid), myByHi(myThid)
195         CALL LEF_ZERO( trair1 ,myThid )            DO bi = myBxLo(myThid), myBxHi(myThid)
196         CALL LEF_ZERO( qrair0 ,myThid )             DO j=1-OLy,sNy+OLy
197         CALL LEF_ZERO( qrair1 ,myThid )              DO i=1-OLx,sNx+OLx
198         CALL LEF_ZERO( solar0 ,myThid )                trair0  (i,j,bi,bj) = 0.
199         CALL LEF_ZERO( solar1 ,myThid )                trair1  (i,j,bi,bj) = 0.
200         CALL LEF_ZERO( uwind0 ,myThid )                qrair0  (i,j,bi,bj) = 0.
201         CALL LEF_ZERO( uwind1 ,myThid )                qrair1  (i,j,bi,bj) = 0.
202         CALL LEF_ZERO( vwind0 ,myThid )                solar0  (i,j,bi,bj) = 0.
203         CALL LEF_ZERO( vwind1 ,myThid )                solar1  (i,j,bi,bj) = 0.
204        ENDIF                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
219               ENDDO
220              ENDDO
221             ENDDO
222            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
   
        _BEGIN_MASTER(myThid)  
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         _END_MASTER(myThid)  C       end of loading new fields block
345  C          ENDIF
        _EXCH_XY_R4(trair0 , myThid )  
        _EXCH_XY_R4(qrair0 , myThid )  
        _EXCH_XY_R4(solar0 , myThid )  
        _EXCH_XY_R4(uwind0 , myThid )  
        _EXCH_XY_R4(vwind0 , myThid )  
        _EXCH_XY_R4(trair1 , myThid )  
        _EXCH_XY_R4(qrair1 , myThid )  
        _EXCH_XY_R4(solar1 , myThid )  
        _EXCH_XY_R4(uwind1 , myThid )  
        _EXCH_XY_R4(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_RS( SolarFile,solar,' ',0,myThid )           CALL READ_FLD_XY_RL( SolarFile,' ',solar,0,myThid )
397        ELSE          ELSE
398        DO bj=1,nSy           DO bj = myByLo(myThid), myByHi(myThid)
399           DO bi=1,nSx            DO bi = myBxLo(myThid), myBxHi(myThid)
400              DO j=1,sNy             DO j=1,sNy
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
                ENDDO  
405              ENDDO              ENDDO
406               ENDDO
407              ENDDO
408           ENDDO           ENDDO
409        ENDDO          ENDIF
410        _EXCH_XY_RS(solar, mythid)          _EXCH_XY_RL(solar, myThid)
       ENDIF  
411          IF ( TrFile .NE. ' '  ) THEN          IF ( TrFile .NE. ' '  ) THEN
412           CALL READ_FLD_XY_RS( TrFile,tr,' ',0,myThid )           CALL READ_FLD_XY_RL( TrFile,' ',tr,0,myThid )
413        ELSE          ELSE
414        DO bj=1,nSy           DO bj = myByLo(myThid), myByHi(myThid)
415           DO bi=1,nSx            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              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  
                ENDDO  
422              ENDDO              ENDDO
423               ENDDO
424              ENDDO
425           ENDDO           ENDDO
       ENDDO  
       _EXCH_XY_RS(TR, mythid)  
426          ENDIF          ENDIF
427  c do specific humidity          _EXCH_XY_RL(Tr, myThid)
428          IF ( QrFile .NE. ' '  ) THEN  
429           CALL READ_FLD_XY_RS( QrFile,qr,' ',0,myThid )  C do specific humidity
430        ELSE          IF ( QrFile .NE. ' '.AND..NOT. useRelativeHumidity  ) THEN
431  c default specific humidity profile to 80% relative humidity           CALL READ_FLD_XY_RL( QrFile,' ',qr,0,myThid )
432        DO bj=1,nSy          ELSEIF ( QrrelFile .NE. ' '.AND.useRelativeHumidity) THEN
433           DO bi=1,nSx           CALL READ_FLD_XY_RL( QrrelFile,' ',rref,0,myThid )
434              DO j=1,sNy           CALL CHEAPAML_CONVERT_HUM(rref, Tr, qr,myThid )
435                 DO i=1,sNx  
436            ELSE
437    C default specific humidity profile to 80% relative humidity
438             DO bj = myByLo(myThid), myByHi(myThid)
439              DO bi = myBxLo(myThid), myBxHi(myThid)
440               DO j=1,sNy
441                DO i=1,sNx
442  c                  jG = myYGlobalLo-1+(bj-1)*sNy+j  c                  jG = myYGlobalLo-1+(bj-1)*sNy+j
443                    local = Tr(i,j,bi,bj)+273.16d0                    local = Tr(i,j,bi,bj)+273.16d0
444                ssqa = ssq0*exp( lath*(ssq1-ssq2/local)) / p0                ssqa = ssq0*exp( lath*(ssq1-ssq2/local)) / p0
445                    qr(i,j,bi,bj) = 0.8d0*ssqa                    qr(i,j,bi,bj) = 0.8d0*ssqa
                ENDDO  
446              ENDDO              ENDDO
447               ENDDO
448              ENDDO
449           ENDDO           ENDDO
       ENDDO  
       _EXCH_XY_RS(qr, mythid)  
450          ENDIF          ENDIF
451            _EXCH_XY_RL(qr, myThid)
452          IF ( UWindFile .NE. ' '  ) THEN          IF ( UWindFile .NE. ' '  ) THEN
453           CALL READ_FLD_XY_RS( UWindFile,uwind,' ',0,myThid )           CALL READ_FLD_XY_RL( UWindFile,' ',uwind,0,myThid )
454        ELSE          ELSE
455        DO bj=1,nSy           DO bj = myByLo(myThid), myByHi(myThid)
456           DO bi=1,nSx            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
                ENDDO  
466              ENDDO              ENDDO
467               ENDDO
468              ENDDO
469           ENDDO           ENDDO
       ENDDO  
       _EXCH_XY_RS(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_RS( VWindFile,vwind,' ',0,myThid )           CALL READ_FLD_XY_RL( VWindFile,' ',vwind,0,myThid )
474        ELSE          ELSE
475        DO bj=1,nSy           DO bj = myByLo(myThid), myByHi(myThid)
476           DO bi=1,nSx            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
                ENDDO  
481              ENDDO              ENDDO
482               ENDDO
483              ENDDO
484           ENDDO           ENDDO
       ENDDO  
       _EXCH_XY_RS(vwind, mythid)  
485          ENDIF          ENDIF
486            _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.2  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22