/[MITgcm]/MITgcm/pkg/aim/aim_do_atmos_physics.F
ViewVC logotype

Diff of /MITgcm/pkg/aim/aim_do_atmos_physics.F

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

revision 1.1 by cnh, Fri Jan 26 00:14:31 2001 UTC revision 1.2 by adcroft, Fri Feb 2 21:36:29 2001 UTC
# Line 0  Line 1 
1    C $Header$
2    C     $Name$
3    
4    #include "AIM_OPTIONS.h"
5    
6          SUBROUTINE AIM_DO_ATMOS_PHYSICS( phi_hyd, currentTime, myThid )
7    C     /==================================================================\
8    C     | S/R AIM_DO_ATMOS_PHYSICS                                         |
9    C     |==================================================================|
10    C     | Interface interface between atmospheric physics package and the  |
11    C     | dynamical model.                                                 |
12    C     | Routine calls physics pacakge after mapping model variables to   |
13    C     | the package grid. Package should derive and set tendency terms   |
14    C     | which can be included as external forcing terms in the dynamical |
15    C     | tendency routines. Packages should communicate this information  |
16    C     | through common blocks.                                           |
17    C     \==================================================================/
18    
19    C     -------------- Global variables ------------------------------------
20    C     Physics package
21    #include "atparam.h"
22    #include "atparam1.h"
23          INTEGER NGP
24          INTEGER NLON
25          INTEGER NLAT
26          INTEGER NLEV
27          PARAMETER ( NLON=IX, NLAT=IL, NLEV=KX, NGP=NLON*NLAT )
28    #include "com_physvar.h"
29    #include "com_forcing1.h"
30    #include "Lev_def.h"
31    C     MITgcm
32    #include "EEPARAMS.h"
33    #include "PARAMS.h"
34    #include "DYNVARS.h"
35    #include "CG2D.h"
36    #include "GRID.h"
37    
38    C     -------------- Routine arguments -----------------------------------
39          _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
40          _RL currentTime
41    
42    #ifdef ALLOW_AIM
43    C     -------------- Local variables -------------------------------------
44    C     I,J,K,I2,J2   - Loop counters
45    C     tYear         - Fraction into year
46    C     mnthIndex     - Current month
47    C     prevMnthIndex - Month last time this routine was called.
48    C     tmp4          - I/O buffer ( 32-bit precision )
49    C     fNam          - Work space for file names
50    C     mnthNam       - Month strings
51    C     hInital       - Initial height of pressure surfaces (m)
52    C     pSurfs        - Pressure surfaces (Pa)
53    C     Katm          - Atmospheric K index
54          INTEGER I
55          INTEGER I2
56          INTEGER J
57          INTEGER J2
58          INTEGER K
59          INTEGER IG0
60          INTEGER JG0
61          REAL    tYear
62          INTEGER mnthIndex
63          INTEGER prevMnthIndex
64          DATA    prevMnthIndex / 0 /
65          SAVE    prevMnthIndex
66          LOGICAL FirstCall
67          DATA    FirstCall /.TRUE./
68          SAVE    FirstCall
69          LOGICAL CALLFirst
70          DATA    CALLFirst /.TRUE./
71          SAVE    CALLFirst
72          INTEGER nxIo
73          INTEGER nyIo
74          PARAMETER ( nxIo = 128, nyIo = 64 )
75          Real*4  tmp4(nxIo,nyIo)
76          CHARACTER*16 fNam
77          CHARACTER*3 mnthNam(12)
78          DATA mnthNam /
79         & 'jan', 'feb', 'mar', 'apr', 'may', 'jun',
80         & 'jul', 'aug', 'sep', 'oct', 'nov', 'dec' /
81          SAVE mnthNam
82          REAL    hInitial(Nr)
83          REAL    hInitialW(Nr)
84          DATA    hInitial / 418.038,2038.54,5296.88,10090.02,17338.0/
85          SAVE    hInitial
86          DATA    hInitialW / 0., 1657.54, 4087.75, 8050.96,15090.4 /
87          REAL    pSurfs(Nr)
88          DATA    pSurfs   / 950.D2,775.D2, 500.D2,  250.D2, 75.D2 /
89          SAVE    pSurfs
90          REAL    pSurfw(Nr)
91          DATA    pSurfw   /1000.D2, 900.D2, 650.D2, 350.D2, 150.D2  /
92          SAVE    pSurfw
93          REAL    RD
94          REAL    CPAIR
95          REAL    RhoG1(sNx*sNy,Nr)
96          INTEGER npasdt
97          DATA npasdt /0/
98          SAVE npasdt
99          REAL Soilqmax
100          REAL phiTotal(sNx,sNy,Nr)
101          _RL  phiTCount
102          _RL  phiTSum
103          _RL  ans
104          real pvoltotNiv5
105          SAVE pvoltotNiv5
106          real ptotalNiv5
107          INTEGER bi, bj
108          INTEGER Katm
109    C
110          pGround = 1.D5
111          CPAIR   = 1004
112          RD      =  287
113    
114    C     Assume only one tile per proc. for now
115          bi  = 1
116          bj  = 1
117          IG0 = myXGlobalLo
118          JG0 = myYGlobalLo
119    
120    C      
121    C     Physics package works with sub-domains 1:sNx,1:sNy,1:Nr.
122    C     Internal index mapping is linear in X and Y with a second
123    C     dimension for the vertical.
124    
125    C     Adjustment for heave due to mean heating/cooling
126    C     ( I don't think the old formula was strictly "correct" for orography
127    C       but I have implemented it as was for now. As implemented
128    C       the mean heave of the bottom (K=Nr) level is calculated rather than
129    C       the mean heave of the base of the atmosphere. )
130          phiTCount = 0.
131          phiTSum   = 0.
132          DO K=1,Nr
133          DO J=1,sNy
134           DO I=1,sNx
135            phiTotal(I,J,K)  = cg2d_x(i,j,bi,bj)
136            phiTCount        = phiTCount + hFacC(i,j,Nr,bi,bj)
137           ENDDO
138          ENDDO
139          ENDDO
140          DO K=1,Nr
141           DO J=1,sNy
142            DO I=1,sNx
143             phiTotal(I,J,K) = phiTotal(I,J,K) +
144         &   recip_rhoConst*(phi_hyd(i,j,k,bi,bj))
145            ENDDO
146           ENDDO
147          ENDDO
148          DO J=1,sNy
149           DO I=1,sNx
150            phiTSum = phiTSum + phiTotal(I,J,Nr)
151           ENDDO
152          ENDDO
153          ans = phiTCount
154    C     _GLOBAL_SUM_R8( phiTCount, myThid )
155          phiTcount = ans
156          ans = phiTSum
157    C     _GLOBAL_SUM_R8( phiTSum, myThid )
158          phiTSum = ans
159    C     ptotalniv5=phiTSum/phiTCount
160          ptotalniv5=0.
161    
162    C     Note the mapping here is only valid for one tile
163    C     per proc.
164          DO K = 1, Nr
165           DO J = 1, sNy
166            DO I = 1, sNx
167             I2 = (sNx)*(J-1)+I
168             Katm = _KD2KA( K )
169             UG1(I2,Katm)   = 0.5*(uVel(I,J,K,bi,bj)+uVel(I+1,J,K,bi,bj))
170             VG1(I2,Katm)   = 0.5*(vVel(I,J,K,bi,bj)+vVel(I,J+1,K,bi,bj))
171    C        Phyiscs works with temperature - not potential temp.
172             TG1(I2,Katm)   = theta(I,J,K,bi,bj)/((pGround/pSurfs(K))**(RD/CPAIR))
173             QG1(I2,Katm)   = salt(I,J,K,bi,bj)
174             PHIG1(I2,Katm) = (phiTotal(I,J,K)- ptotalniv5 ) + gravity*Hinitial(k)
175             if(hFacC(i,j,k,bi,bj).eq.1.) then
176               RHOG1(I2,Katm) = pSurfs(K)/RD/TG1(I2,Katm)
177             else
178               RHOG1(I2,Katm)=0.
179             endif
180            ENDDO
181           ENDDO
182          ENDDO
183    
184    C      
185    C     Set geopotential surfaces
186    C     -------------------------
187          DO J=1,sNy
188           DO I=1,sNx
189            I2 = (sNx)*(J-1)+I
190            IF ( Nlevxy(I2) .NE. 0 ) THEN
191             PHI0(I2) = gravity*Hinitialw(Nlevxy(I2))
192            ELSE
193             PHI0(I2) = 0.
194            ENDIF
195           ENDDO
196          ENDDO
197    C
198    C     Physics package works with log of surface pressure
199    C     Get surface pressure from pbot-dpref/dz*Z'
200          DO J=1,sNy
201           DO I=1,sNx
202            I2 = (sNx)*(J-1)+I
203            IF ( Nlevxy(I2) .NE. 0 ) THEN
204             PNLEVW(I2) = PsurfW(Nlevxy(I2))/pGround
205            ELSE
206    C        Dummy value for land
207             PNLEVW(I2) = PsurfW(1)/pGround
208            ENDIF
209            PSLG1(I2) = 0.
210           ENDDO
211          ENDDO
212    cch      write(0,*)  '(PNLEVW(I2),I2=257,384)'
213    cch      write(0,*)  (PNLEVW(I2),I2=257,384)
214    C
215    C
216    C     Physics package needs to know time of year as a fraction
217          tYear = currentTime/(86400.*360.) -
218         &        FLOAT(INT(currentTime/(86400.*360.)))
219    C
220    C     Load external data needed by physics package
221    C     1. Albedo
222    C     2. Soil moisture
223    C     3. Surface temperatures
224    C     4. Snow depth            - assume no snow for now
225    C     5. Sea ice               - assume no sea ice for now
226    C     6. Land sea mask         - infer from exact zeros in soil moisture dataset
227    C     7. Surface geopotential  - to be done when orography is in
228    C                                dynamical kernel. Assume 0. for now.
229          mnthIndex = INT(tYear*12.)+1
230          IF ( mnthIndex .NE. prevMnthIndex .OR.
231         &     FirstCall ) THEN
232           prevMnthIndex = mnthIndex
233    C      Read in surface albedo data (input is in % 0-100 )
234    C      scale to give fraction between 0-1 for Francos package.
235    CequChan       WRITE(fNam,'(A,A,A)' ) 'salb.',mnthNam(mnthIndex),'.sun.b'
236    CequChan       OPEN(1,FILE=fNam(1:14),STATUS='old',FORM='unformatted')
237    CequChan       READ(1) tmp4
238    CequChan       CLOSE(1)
239    CequChan       DO J=1,nYio
240    CequChan        DO I=1,nXio
241    CequChan         tmp4(I,J) = tmp4(I,J)/100.
242    CequChan        ENDDO
243    CequChan       ENDDO
244           DO J=1,sNy
245            DO I=1,sNx
246             I2 = (sNx)*(J-1)+I
247             alb0(I2) = 0.
248    CequChan         IF ( IG0+I-1 .LE. nxIo .AND. JG0+J-1 .LE. nyIo ) THEN
249    CequChan          alb0(I2) = tmp4(IG0+I-1,JG0+J-1)
250    CequChan         ENDIF
251            ENDDO
252           ENDDO
253    C      Read in surface temperature data (input is in absolute temperature)
254    CequChan       WRITE(fNam,'(A,A,A)' ) 'tsurf.',mnthNam(mnthIndex),'.sun.b'
255    CequChan       OPEN(1,FILE=fNam(1:15),STATUS='old',FORM='unformatted')
256    CequChan       READ(1) tmp4
257    CequChan       CLOSE(1)
258           DO J=1,sNy
259            DO I=1,sNx
260             I2 = (sNx)*(J-1)+I
261             sst1(I2) = 300.
262             stl1(I2) = 300.
263    CequChan         IF ( IG0+I-1 .LE. nxIo .AND. JG0+J-1 .LE. nyIo ) THEN
264    CequChan          sst1(I2) = tmp4(IG0+I-1,JG0+J-1)
265    CequChan          stl1(I2) = tmp4(IG0+I-1,JG0+J-1)
266    CequChan         ENDIF
267    caja     IF ( I .GE. 64-10 .AND. I .LE. 65+10 ) THEN
268    caja      sst1(I2) = 310.
269    caja      stl1(I2) = 310.
270    caja     ENDIF
271    caja     IF ( I .GE. 64-10 .AND. I .LE. 65+10 ) THEN
272    caja      sst1(I2) = 300.+10.*exp( -((float(I)-64.5)/5.)**2 )
273    caja      stl1(I2) = sst1(I2)
274    caja     ENDIF
275    c_jmc: should not be part of the AIM package :
276             sst1(I2) = 300.+10.*exp( -((float(I)-64.5)/25.)**2 )
277             stl1(I2) = sst1(I2)
278            ENDDO
279           ENDDO
280    C
281    C      Read in soil moisture data (input is in cm in bucket of depth 20cm. )
282    C??? NOT CLEAR  scale for bucket depth of 75mm which is what Franco uses.
283    CequChan       WRITE(fNam,'(A,A,A)' ) 'smoist.',mnthNam(mnthIndex),'.sun.b'
284    CequChan       OPEN(1,FILE=fNam(1:16),STATUS='old',FORM='unformatted')
285    CequChan       READ(1) tmp4
286    CequChan       CLOSE(1)
287    CequChan       WRITE(0,*) ' Read file ', fNam(1:16), IG0, JG0
288    cdj       tmp4 = (tmp4*7.5/20.)*10.
289           DO J=1,sNy
290            DO I=1,sNx
291             I2 = (sNx)*(J-1)+I
292             soilq1(I2) = 0.
293    CequChan         IF ( IG0+I-1 .LE. nxIo .AND. JG0+J-1 .LE. nyIo ) THEN
294    CequChan          soilq1(I2) = tmp4(IG0+I-1,JG0+J-1)
295    CequChan         ENDIF
296            ENDDO
297           ENDDO
298    cdj      Soilqmax=MAxval(soilq1)
299           Soilqmax=20.
300    cdj      if(Soilqmax.ne.0.) then
301           DO J=1,sNy
302            DO I=1,sNx
303             I2 = (sNx)*(J-1)+I
304    CequChan         soilq1(I2)=soilq1(I2)/Soilqmax
305             soilq1(I2) = 1.
306            ENDDO
307           ENDDO
308    cdj      endif
309          ENDIF
310    C
311          IF ( FirstCall ) THEN
312    C      Set snow depth, sea ice to zero for now
313    C      Land-sea mask ( figure this out from where soil moisture is exactly zero ).
314           DO J=1,sNy
315            DO I=1,sNx
316             I2 = (sNx)*(J-1)+I
317             fMask1(I2) = 1.
318             IF ( soilq1(I2) .EQ. 0. ) fMask1(I2) = 0.
319             oice1(I2) = 0.
320             snow1(I2) = 0.
321            ENDDO
322           ENDDO
323    C      open(77,file='lsmask',form='unformatted')
324    C      write(77) fmask1
325    C      close(77)
326          ENDIF
327    C
328    C Addition may 15 . Reset humidity to 0. if negative
329    C ---------------------------------------------------
330          DO K=1,Nr
331           DO J=1-OLy,sNy+OLy
332            DO I=1-Olx,sNx+OLx
333             IF ( salt(i,j,k,bi,bj) .LT. 0. .OR. K .EQ. Nr ) THEN
334              salt(i,j,k,bi,bj) = 0.
335             ENDIF
336            ENDDO
337           ENDDO
338          ENDDO
339    C
340          CALL PDRIVER( tYear )
341    
342    #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE
343    C     Calculate diagnostics for AIM
344          CALL AIM_CALC_DIAGS( bi, bj, currentTime, myThid )
345    #endif /* INCLUDE_DIAGNOSTICS_INTERFACE_CODE */
346    C
347          FirstCall = .FALSE.
348    C
349    #endif /* ALLOW_AIM */
350    
351          RETURN
352          END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22