/[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.9 by jmc, Tue Jan 8 22:33:10 2002 UTC revision 1.10 by jmc, Fri Sep 27 20:05:11 2002 UTC
# Line 17  C     | which can be included as externa Line 17  C     | which can be included as externa
17  C     | tendency routines. Packages should communicate this information  |  C     | tendency routines. Packages should communicate this information  |
18  C     | through common blocks.                                           |  C     | through common blocks.                                           |
19  C     \==================================================================/  C     \==================================================================/
20        IMPLICIT rEAL*8 (A-H,O-Z)        IMPLICIT NONE
21    
22  C     -------------- Global variables ------------------------------------  C     -------------- Global variables ------------------------------------
23  C     Physics package  C-- size for MITgcm & Physics package :
24  #include "atparam.h"  #include "AIM_SIZE.h"
 #include "atparam1.h"  
       INTEGER NGP  
       INTEGER NLON  
       INTEGER NLAT  
       INTEGER NLEV  
       PARAMETER ( NLON=IX, NLAT=IL, NLEV=KX, NGP=NLON*NLAT )  
25    
26  C     MITgcm  C-- MITgcm
27  #include "EEPARAMS.h"  #include "EEPARAMS.h"
28  #include "PARAMS.h"  #include "PARAMS.h"
29  #include "DYNVARS.h"  #include "DYNVARS.h"
30  #include "GRID.h"  #include "GRID.h"
31  #include "SURFACE.h"  #include "SURFACE.h"
 #include "AIM_FFIELDS.h"  
32    
33  C     Physics package  C-- Physics package
34    #include "AIM_FFIELDS.h"
35    #include "AIM_GRID.h"
36  #include "com_physvar.h"  #include "com_physvar.h"
37  #include "com_forcing1.h"  #include "com_forcing1.h"
 #include "Lev_def.h"  
38    
39  C     -------------- Routine arguments -----------------------------------  C     -------------- Routine arguments -----------------------------------
40        _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 50  C     -------------- Routine arguments - Line 44  C     -------------- Routine arguments -
44    
45  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
46  C     -------------- Local variables -------------------------------------  C     -------------- Local variables -------------------------------------
47  C     I,J,K,I2,J2   - Loop counters  C     I,J,K,I2      - Loop counters
48  C     tYear         - Fraction into year  C     tYear         - Fraction into year
 C     mnthIndex     - Current month  
 C     prevMnthIndex - Month last time this routine was called.  
 C     tmp4          - I/O buffer ( 32-bit precision )  
 C     fNam          - Work space for file names  
 C     mnthNam       - Month strings  
49  C     hInital       - Initial height of pressure surfaces (m)  C     hInital       - Initial height of pressure surfaces (m)
50  C     pSurfs        - Pressure surfaces (Pa)  C     pSurfs        - Pressure surfaces (Pa)
51  C     Katm          - Atmospheric K index  C     Katm          - Atmospheric K index
52        INTEGER I        INTEGER I
53        INTEGER I2        INTEGER I2
54        INTEGER J        INTEGER J
       INTEGER J2  
55        INTEGER K        INTEGER K
56        INTEGER IG0        _RL     tYear
57        INTEGER JG0        _RL     hInitial(Nr)
58        REAL    tYear        _RL     hInitialW(Nr)
       INTEGER mnthIndex  
       INTEGER prevMnthIndex  
       DATA    prevMnthIndex / 0 /  
       SAVE    prevMnthIndex  
       LOGICAL FirstCall  
       DATA    FirstCall /.TRUE./  
       SAVE    FirstCall  
       LOGICAL CALLFirst  
       DATA    CALLFirst /.TRUE./  
       SAVE    CALLFirst  
       INTEGER nxIo  
       INTEGER nyIo  
       PARAMETER ( nxIo = 128, nyIo = 64 )  
       Real*4  tmp4(nxIo,nyIo)  
       CHARACTER*16 fNam  
       CHARACTER*3 mnthNam(12)  
       DATA mnthNam /  
      & 'jan', 'feb', 'mar', 'apr', 'may', 'jun',  
      & 'jul', 'aug', 'sep', 'oct', 'nov', 'dec' /  
       SAVE mnthNam  
       REAL    hInitial(Nr)  
       REAL    hInitialW(Nr)  
59        DATA    hInitial / 17338.0,10090.02,5296.88,2038.54,418.038/        DATA    hInitial / 17338.0,10090.02,5296.88,2038.54,418.038/
60        SAVE    hInitial        SAVE    hInitial
61        DATA    hInitialW / 15090.4, 8050.96, 4087.75, 1657.54, 0. /        DATA    hInitialW / 15090.4, 8050.96, 4087.75, 1657.54, 0. /
62        REAL    pSurfs(Nr)        _RL     pSurfs(Nr)
63        DATA    pSurfs   / 75.D2, 250.D2, 500.D2, 775.D2, 950.D2 /        DATA    pSurfs   / 75.D2, 250.D2, 500.D2, 775.D2, 950.D2 /
64        SAVE    pSurfs        SAVE    pSurfs
65        REAL    pSurfw(Nr)        _RL     pSurfw(Nr)
66        DATA    pSurfw / 150.D2, 350.D2, 650.D2, 900.D2, 1000.D2 /        DATA    pSurfw / 150.D2, 350.D2, 650.D2, 900.D2, 1000.D2 /
67        SAVE    pSurfw        SAVE    pSurfw
68        REAL    RD        _RL     RD
69        REAL    CPAIR        _RL     CPAIR
70        REAL    RhoG1(sNx*sNy,Nr)        _RL     pGround
71        INTEGER npasdt        _RL     RhoG1(sNx*sNy,Nr)
72        DATA npasdt /0/  c     _RL  phiTotal(sNx,sNy,Nr)
73        SAVE npasdt  c     _RL  phiTCount
74        REAL Soilqmax  c     _RL  phiTSum
75        REAL phiTotal(sNx,sNy,Nr)  c     _RL  ans
76        _RL  phiTCount  c     _RL  pvoltotNiv5
77        _RL  phiTSum  c     SAVE pvoltotNiv5
78        _RL  ans  c     _RL  ptotalNiv5
       real pvoltotNiv5  
       SAVE pvoltotNiv5  
       real ptotalNiv5  
79        INTEGER Katm        INTEGER Katm
80    
 C  
81        pGround = 1.D5        pGround = 1.D5
82        CPAIR   = 1004        CPAIR   = 1004
83        RD      =  287        RD      =  287
84    
85        CALL AIM_DYN2AIM( bi, bj, currentTime, myThid )        CALL AIM_DYN2AIM( bi, bj, currentTime, myThid )
86    
 C     Assume only one tile per proc. for now  
       IG0 = myXGlobalLo+(bi-1)*sNx  
       JG0 = myYGlobalLo+(bj-1)*sNy  
   
 C        
87  C     Physics package works with sub-domains 1:sNx,1:sNy,1:Nr.  C     Physics package works with sub-domains 1:sNx,1:sNy,1:Nr.
88  C     Internal index mapping is linear in X and Y with a second  C     Internal index mapping is linear in X and Y with a second
89  C     dimension for the vertical.  C     dimension for the vertical.
# Line 136  C     ( I don't think the old formula wa Line 93  C     ( I don't think the old formula wa
93  C       but I have implemented it as was for now. As implemented  C       but I have implemented it as was for now. As implemented
94  C       the mean heave of the bottom (K=Nr) level is calculated rather than  C       the mean heave of the bottom (K=Nr) level is calculated rather than
95  C       the mean heave of the base of the atmosphere. )  C       the mean heave of the base of the atmosphere. )
96        phiTCount = 0.  c     phiTCount = 0.
97        phiTSum   = 0.  c     phiTSum   = 0.
98        DO K=1,Nr  c     DO K=1,Nr
99        DO J=1,sNy  c     DO J=1,sNy
100         DO I=1,sNx  c      DO I=1,sNx
101          phiTotal(I,J,K)  = etaN(i,j,bi,bj)  c       phiTotal(I,J,K)  = etaN(i,j,bi,bj)
102          phiTCount        = phiTCount + hFacC(i,j,Nr,bi,bj)  c       phiTCount        = phiTCount + hFacC(i,j,Nr,bi,bj)
103         ENDDO  c      ENDDO
104        ENDDO  c     ENDDO
105        ENDDO  c     ENDDO
106        DO K=1,Nr  c     DO K=1,Nr
107         DO J=1,sNy  c      DO J=1,sNy
108          DO I=1,sNx  c       DO I=1,sNx
109           phiTotal(I,J,K) = phiTotal(I,J,K) +  c        phiTotal(I,J,K) = phiTotal(I,J,K) +
110       &   recip_rhoConst*(phi_hyd(i,j,k))  c    &   recip_rhoConst*(phi_hyd(i,j,k))
111          ENDDO  c       ENDDO
112         ENDDO  c      ENDDO
113        ENDDO  c     ENDDO
114        DO J=1,sNy  c     DO J=1,sNy
115         DO I=1,sNx  c      DO I=1,sNx
116          phiTSum = phiTSum + phiTotal(I,J,Nr)  c       phiTSum = phiTSum + phiTotal(I,J,Nr)
117         ENDDO  c      ENDDO
118        ENDDO  c     ENDDO
119        ans = phiTCount  c     ans = phiTCount
120  C     _GLOBAL_SUM_R8( phiTCount, myThid )  C     _GLOBAL_SUM_R8( phiTCount, myThid )
121        phiTcount = ans  c     phiTcount = ans
122        ans = phiTSum  c     ans = phiTSum
123  C     _GLOBAL_SUM_R8( phiTSum, myThid )  C     _GLOBAL_SUM_R8( phiTSum, myThid )
124        phiTSum = ans  c     phiTSum = ans
125  C     ptotalniv5=phiTSum/phiTCount  C     ptotalniv5=phiTSum/phiTCount
126        ptotalniv5=0.  c     ptotalniv5=0.
127    
128  #ifndef OLD_AIM_INTERFACE  #ifndef OLD_AIM_INTERFACE
129  c_jmc: Because AIM physics LSC is not applied in the stratosphere (top level),  C_jmc: Because AIM physics LSC is not applied in the stratosphere (top level),
130  c      ==> move water wapor from the stratos to the surface level.  C      ==> move water wapor from the stratos to the surface level.
131        DO J = 1-Oly, sNy+Oly        DO J = 1-Oly, sNy+Oly
132         DO I = 1-Olx, sNx+Olx         DO I = 1-Olx, sNx+Olx
133          k = ksurfC(i,j,bi,bj)          k = ksurfC(i,j,bi,bj)
# Line 199  C        Physics works with temperature Line 156  C        Physics works with temperature
156  #else  #else
157           QG1(I2,Katm,myThid)   = MAX(salt(I,J,K,bi,bj), 0. _d 0)           QG1(I2,Katm,myThid)   = MAX(salt(I,J,K,bi,bj), 0. _d 0)
158  #endif  #endif
159           PHIG1(I2,Katm,myThid) = (phiTotal(I,J,K)- ptotalniv5 )           PHIG1(I2,Katm,myThid) = gravity*Hinitial(Katm)
160       &                  + gravity*Hinitial(Katm)  c    &                         +(phiTotal(I,J,K)- ptotalniv5 )
 C *NOTE* Fix me for lopped cells <== done !  
161           IF (maskC(i,j,k,bi,bj).EQ.1.) THEN           IF (maskC(i,j,k,bi,bj).EQ.1.) THEN
162             RHOG1(I2,Katm) = pSurfs(Katm)/RD/TG1(I2,Katm,myThid)             RHOG1(I2,Katm) = pSurfs(Katm)/RD/TG1(I2,Katm,myThid)
163           ELSE           ELSE
# Line 212  C *NOTE* Fix me for lopped cells <== don Line 168  C *NOTE* Fix me for lopped cells <== don
168        ENDDO        ENDDO
169    
170  c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
171  c_jmc: add square of surface wind speed (center of C grid) = 2 * KE_surf  C_jmc: add square of surface wind speed (center of C grid) = 2 * KE_surf
172        DO J = 1, sNy        DO J = 1, sNy
173         DO I = 1, sNx         DO I = 1, sNx
174          I2 = I+(J-1)*sNx          I2 = I+(J-1)*sNx
# Line 273  C        Dummy value for land Line 229  C        Dummy value for land
229          PSLG1(I2,myThid) = 0.          PSLG1(I2,myThid) = 0.
230         ENDDO         ENDDO
231        ENDDO        ENDDO
 cch      write(0,*)  '(PNLEVW(I2),I2=257,384)'  
 cch      write(0,*)  (PNLEVW(I2),I2=257,384)  
232  C  C
233  C  C
234  C     Physics package needs to know time of year as a fraction  C     Physics package needs to know time of year as a fraction
# Line 291  C     5. Sea ice               - assume Line 245  C     5. Sea ice               - assume
245  C     6. Land sea mask         - infer from exact zeros in soil moisture dataset  C     6. Land sea mask         - infer from exact zeros in soil moisture dataset
246  C     7. Surface geopotential  - to be done when orography is in  C     7. Surface geopotential  - to be done when orography is in
247  C                                dynamical kernel. Assume 0. for now.  C                                dynamical kernel. Assume 0. for now.
248        mnthIndex = INT(tYear*12.)+1  
249  C_cnh01      IF ( mnthIndex .NE. prevMnthIndex .OR.  C      Load in surface albedo data (in [0,1]) from aim_albedo to alb0 :
 C_cnh01     &     FirstCall ) THEN  
 C_cnh01       prevMnthIndex = mnthIndex  
 C      Read in surface albedo data (input is in % 0-100 )  
 C      scale to give fraction between 0-1 for Francos package.  
 C      WRITE(fNam,'(A,A,A)' ) 'salb.',mnthNam(mnthIndex),'.sun.b'  
 C      OPEN(1,FILE=fNam(1:14),STATUS='old',FORM='unformatted')  
 C      READ(1) tmp4  
 C      CLOSE(1)  
 C      DO J=1,nYio  
 C       DO I=1,nXio  
 C        tmp4(I,J) = aim_albedo(I,J)/100.  
 C       ENDDO  
 C      ENDDO  
250         DO J=1,sNy         DO J=1,sNy
251          DO I=1,sNx          DO I=1,sNx
252           I2 = (sNx)*(J-1)+I           I2 = (sNx)*(J-1)+I
253           alb0(I2,myThid) = 0.           alb0(I2,myThid) = 0.
 c        alb0(I2,myThid) = aim_albedo(I,J,bi,bj)/100.  
254           alb0(I2,myThid) = aim_albedo(I,J,bi,bj)           alb0(I2,myThid) = aim_albedo(I,J,bi,bj)
255          ENDDO          ENDDO
256         ENDDO         ENDDO
257  C      Read in surface temperature data (input is in absolute temperature)  C      Load in surface temperature data from aim_surfTemp to stl1 & sst1 :
 C      WRITE(fNam,'(A,A,A)' ) 'tsurf.',mnthNam(mnthIndex),'.sun.b'  
 C      OPEN(1,FILE=fNam(1:15),STATUS='old',FORM='unformatted')  
 C      READ(1) tmp4  
 C      CLOSE(1)  
258         DO J=1,sNy         DO J=1,sNy
259          DO I=1,sNx          DO I=1,sNx
260           I2 = (sNx)*(J-1)+I           I2 = (sNx)*(J-1)+I
# Line 329  C      CLOSE(1) Line 265  C      CLOSE(1)
265          ENDDO          ENDDO
266         ENDDO         ENDDO
267  C  C
268  C      Read in soil moisture data (input is in cm in bucket of depth 20cm. )  C      Load in soil moisture data (in [0,1]) from aim_soilMoisture to soilq1 :
269  C??? NOT CLEAR  scale for bucket depth of 75mm which is what Franco uses.  C??? NOT CLEAR  scale for bucket depth of 75mm which is what Franco uses.
 C      WRITE(fNam,'(A,A,A)' ) 'smoist.',mnthNam(mnthIndex),'.sun.b'  
 C      OPEN(1,FILE=fNam(1:16),STATUS='old',FORM='unformatted')  
 C      READ(1) tmp4  
 C      CLOSE(1)  
 C      WRITE(0,*) ' Read file ', fNam(1:16), IG0, JG0  
 cdj       tmp4 = (tmp4*7.5/20.)*10.  
270         DO J=1,sNy         DO J=1,sNy
271          DO I=1,sNx          DO I=1,sNx
272           I2 = (sNx)*(J-1)+I           I2 = (sNx)*(J-1)+I
273           soilq1(I2,myThid) = 0.           soilq1(I2,myThid) = 0.
 c        soilq1(I2,myThid) = aim_soilMoisture(I,J,bi,bj)/20.  
274           soilq1(I2,myThid) = aim_soilMoisture(I,J,bi,bj)           soilq1(I2,myThid) = aim_soilMoisture(I,J,bi,bj)
275          ENDDO          ENDDO
276         ENDDO         ENDDO
 C_cnh01      ENDIF  
 C  
 C_cnh01      IF ( FirstCall ) THEN  
277  C      Set snow depth, sea ice to zero for now  C      Set snow depth, sea ice to zero for now
278  C      Land-sea mask ( figure this out from where  C      Land-sea mask ( figure this out from where
279  C                      soil moisture is exactly zero ).  C                      soil moisture is exactly zero ).
# Line 363  C                      soil moisture is Line 289  C                      soil moisture is
289  C      open(77,file='lsmask',form='unformatted')  C      open(77,file='lsmask',form='unformatted')
290  C      write(77) fmask1  C      write(77) fmask1
291  C      close(77)  C      close(77)
292  C_cnh01      ENDIF  
 C  
293  C Addition may 15 . Reset humidity to 0. if negative  C Addition may 15 . Reset humidity to 0. if negative
294  C ---------------------------------------------------  C ---------------------------------------------------
295  #ifdef OLD_AIM_INTERFACE  #ifdef OLD_AIM_INTERFACE
# Line 385  C -------------------------------------- Line 310  C --------------------------------------
310  C     Calculate diagnostics for AIM  C     Calculate diagnostics for AIM
311        CALL AIM_CALC_DIAGS( bi, bj, currentTime, myThid )        CALL AIM_CALC_DIAGS( bi, bj, currentTime, myThid )
312  #endif /* ALLOW_TIMEAVE */  #endif /* ALLOW_TIMEAVE */
 C  
       FirstCall = .FALSE.  
313    
314        CALL AIM_AIM2DYN( bi, bj, currentTime, myThid )        CALL AIM_AIM2DYN( bi, bj, currentTime, myThid )
315  C  
316  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
317    
318        RETURN        RETURN
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22