/[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.8.2.1 by adcroft, Tue Feb 26 16:04:48 2002 UTC revision 1.14 by cnh, Mon Aug 1 19:34:57 2005 UTC
# Line 3  C     $Name$ Line 3  C     $Name$
3    
4  #include "AIM_OPTIONS.h"  #include "AIM_OPTIONS.h"
5    
6        SUBROUTINE AIM_DO_ATMOS_PHYSICS( phi_hyd,        SUBROUTINE AIM_DO_PHYSICS( bi, bj, currentTime, myIter, myThid )
7       I                                 bi, bj,                                        
      I                                 currentTime, myThid )  
8  C     /==================================================================\  C     /==================================================================\
9  C     | S/R AIM_DO_ATMOS_PHYSICS                                         |  C     | S/R AIM_DO_ATMOS_PHYSICS                                         |
10  C     |==================================================================|  C     |==================================================================|
# Line 17  C     | which can be included as externa Line 16  C     | which can be included as externa
16  C     | tendency routines. Packages should communicate this information  |  C     | tendency routines. Packages should communicate this information  |
17  C     | through common blocks.                                           |  C     | through common blocks.                                           |
18  C     \==================================================================/  C     \==================================================================/
19        IMPLICIT rEAL*8 (A-H,O-Z)        IMPLICIT NONE
20    
21  C     -------------- Global variables ------------------------------------  C     -------------- Global variables ------------------------------------
22  C     Physics package  C-- size for MITgcm & Physics package :
23  #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 )  
24    
25  C     MITgcm  C-- MITgcm
26  #include "EEPARAMS.h"  #include "EEPARAMS.h"
27  #include "PARAMS.h"  #include "PARAMS.h"
28  #include "DYNVARS.h"  #include "DYNVARS.h"
29  #include "GRID.h"  #include "GRID.h"
30  #include "SURFACE.h"  #include "SURFACE.h"
 #include "AIM_FFIELDS.h"  
31    
32  C     Physics package  C-- Physics package
33    #include "AIM_FFIELDS.h"
34    #include "AIM_GRID.h"
35  #include "com_physvar.h"  #include "com_physvar.h"
36  #include "com_forcing1.h"  #include "com_forcing1.h"
 #include "Lev_def.h"  
37    
38  C     -------------- Routine arguments -----------------------------------  C     -------------- Routine arguments -----------------------------------
39        _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  c     _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
40        _RL currentTime        _RL currentTime
41        INTEGER myThid        INTEGER myIter, myThid
42        INTEGER bi, bj        INTEGER bi, bj
43    
44  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
45  C     -------------- Local variables -------------------------------------  C     -------------- Local variables -------------------------------------
46  C     I,J,K,I2,J2   - Loop counters  C     I,J,K,I2      - Loop counters
47  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  
48  C     hInital       - Initial height of pressure surfaces (m)  C     hInital       - Initial height of pressure surfaces (m)
49  C     pSurfs        - Pressure surfaces (Pa)  C     pSurfs        - Pressure surfaces (Pa)
50  C     Katm          - Atmospheric K index  C     Katm          - Atmospheric K index
51        INTEGER I        INTEGER I
52        INTEGER I2        INTEGER I2
53        INTEGER J        INTEGER J
       INTEGER J2  
54        INTEGER K        INTEGER K
55        INTEGER IG0        _RL     tYear
56        INTEGER JG0        _RL     hInitial(Nr)
57        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)  
58        DATA    hInitial / 17338.0,10090.02,5296.88,2038.54,418.038/        DATA    hInitial / 17338.0,10090.02,5296.88,2038.54,418.038/
59        SAVE    hInitial        SAVE    hInitial
60        DATA    hInitialW / 15090.4, 8050.96, 4087.75, 1657.54, 0. /        DATA    hInitialW / 15090.4, 8050.96, 4087.75, 1657.54, 0. /
61        REAL    pSurfs(Nr)        _RL     pSurfs(Nr)
62        DATA    pSurfs   / 75.D2, 250.D2, 500.D2, 775.D2, 950.D2 /        DATA    pSurfs   / 75.D2, 250.D2, 500.D2, 775.D2, 950.D2 /
63        SAVE    pSurfs        SAVE    pSurfs
64        REAL    pSurfw(Nr)        _RL     pSurfw(Nr)
65        DATA    pSurfw / 150.D2, 350.D2, 650.D2, 900.D2, 1000.D2 /        DATA    pSurfw / 150.D2, 350.D2, 650.D2, 900.D2, 1000.D2 /
66        SAVE    pSurfw        SAVE    pSurfw
67        REAL    RD        _RL     RD
68        REAL    CPAIR        _RL     CPAIR
69        REAL    RhoG1(sNx*sNy,Nr)        _RL     pGround
70        INTEGER npasdt        _RL     RhoG1(sNx*sNy,Nr)
71        DATA npasdt /0/  c     _RL  phiTotal(sNx,sNy,Nr)
72        SAVE npasdt  c     _RL  phiTCount
73        REAL Soilqmax  c     _RL  phiTSum
74        REAL phiTotal(sNx,sNy,Nr)  c     _RL  ans
75        _RL  phiTCount  c     _RL  pvoltotNiv5
76        _RL  phiTSum  c     SAVE pvoltotNiv5
77        _RL  ans  c     _RL  ptotalNiv5
       real pvoltotNiv5  
       SAVE pvoltotNiv5  
       real ptotalNiv5  
78        INTEGER Katm        INTEGER Katm
79    
 C  
80        pGround = 1.D5        pGround = 1.D5
81        CPAIR   = 1004        CPAIR   = 1004
82        RD      =  287        RD      =  287
83    
84        CALL AIM_DYN2AIM( bi, bj, currentTime, myThid )        CALL AIM_DYN2AIM( bi, bj, currentTime, myThid )
85    
 C     Assume only one tile per proc. for now  
       IG0 = myXGlobalLo+(bi-1)*sNx  
       JG0 = myYGlobalLo+(bj-1)*sNy  
   
 C        
86  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.
87  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
88  C     dimension for the vertical.  C     dimension for the vertical.
89    
90  C     Adjustment for heave due to mean heating/cooling  C     Adjustment for heave due to mean heating/cooling
91  C     ( I don't think the old formula was strictly "correct" for orography  C     ( I do not think the old formula was strictly "correct" for orography
92  C       but I have implemented it as was for now. As implemented  C       but I have implemented it as was for now. As implemented
93  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
94  C       the mean heave of the base of the atmosphere. )  C       the mean heave of the base of the atmosphere. )
95        phiTCount = 0.  c     phiTCount = 0.
96        phiTSum   = 0.  c     phiTSum   = 0.
97        DO K=1,Nr  c     DO K=1,Nr
98        DO J=1,sNy  c     DO J=1,sNy
99         DO I=1,sNx  c      DO I=1,sNx
100          phiTotal(I,J,K)  = etaN(i,j,bi,bj)  c       phiTotal(I,J,K)  = etaN(i,j,bi,bj)
101          phiTCount        = phiTCount + hFacC(i,j,Nr,bi,bj)  c       phiTCount        = phiTCount + hFacC(i,j,Nr,bi,bj)
102         ENDDO  c      ENDDO
103        ENDDO  c     ENDDO
104        ENDDO  c     ENDDO
105        DO K=1,Nr  c     DO K=1,Nr
106         DO J=1,sNy  c      DO J=1,sNy
107          DO I=1,sNx  c       DO I=1,sNx
108           phiTotal(I,J,K) = phiTotal(I,J,K) +  c        phiTotal(I,J,K) = phiTotal(I,J,K) +
109       &   recip_rhoConst*(phi_hyd(i,j,k))  c    &   recip_rhoConst*(phi_hyd(i,j,k))
110          ENDDO  c       ENDDO
111         ENDDO  c      ENDDO
112        ENDDO  c     ENDDO
113        DO J=1,sNy  c     DO J=1,sNy
114         DO I=1,sNx  c      DO I=1,sNx
115          phiTSum = phiTSum + phiTotal(I,J,Nr)  c       phiTSum = phiTSum + phiTotal(I,J,Nr)
116         ENDDO  c      ENDDO
117        ENDDO  c     ENDDO
118        ans = phiTCount  c     ans = phiTCount
119  C     _GLOBAL_SUM_R8( phiTCount, myThid )  C     _GLOBAL_SUM_R8( phiTCount, myThid )
120        phiTcount = ans  c     phiTcount = ans
121        ans = phiTSum  c     ans = phiTSum
122  C     _GLOBAL_SUM_R8( phiTSum, myThid )  C     _GLOBAL_SUM_R8( phiTSum, myThid )
123        phiTSum = ans  c     phiTSum = ans
124  C     ptotalniv5=phiTSum/phiTCount  C     ptotalniv5=phiTSum/phiTCount
125        ptotalniv5=0.  c     ptotalniv5=0.
126    
127  #ifndef OLD_AIM_INTERFACE  #ifndef OLD_AIM_INTERFACE
128  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),
129  c      ==> move water wapor from the stratos to the surface level.  C      ==> move water wapor from the stratos to the surface level.
130        DO J = 1-Oly, sNy+Oly        DO J = 1-Oly, sNy+Oly
131         DO I = 1-Olx, sNx+Olx         DO I = 1-Olx, sNx+Olx
132          k = ksurfC(i,j,bi,bj)          k = ksurfC(i,j,bi,bj)
# Line 199  C        Physics works with temperature Line 155  C        Physics works with temperature
155  #else  #else
156           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)
157  #endif  #endif
158           PHIG1(I2,Katm,myThid) = (phiTotal(I,J,K)- ptotalniv5 )           PHIG1(I2,Katm,myThid) = gravity*Hinitial(Katm)
159       &                  + gravity*Hinitial(Katm)  c    &                         +(phiTotal(I,J,K)- ptotalniv5 )
 C *NOTE* Fix me for lopped cells <== done !  
160           IF (maskC(i,j,k,bi,bj).EQ.1.) THEN           IF (maskC(i,j,k,bi,bj).EQ.1.) THEN
161             RHOG1(I2,Katm) = pSurfs(Katm)/RD/TG1(I2,Katm,myThid)             RHOG1(I2,Katm) = pSurfs(Katm)/RD/TG1(I2,Katm,myThid)
162           ELSE           ELSE
# Line 212  C *NOTE* Fix me for lopped cells <== don Line 167  C *NOTE* Fix me for lopped cells <== don
167        ENDDO        ENDDO
168    
169  c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
170  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
171        DO J = 1, sNy        DO J = 1, sNy
172         DO I = 1, sNx         DO I = 1, sNx
173          I2 = I+(J-1)*sNx          I2 = I+(J-1)*sNx
# Line 273  C        Dummy value for land Line 228  C        Dummy value for land
228          PSLG1(I2,myThid) = 0.          PSLG1(I2,myThid) = 0.
229         ENDDO         ENDDO
230        ENDDO        ENDDO
 cch      write(0,*)  '(PNLEVW(I2),I2=257,384)'  
 cch      write(0,*)  (PNLEVW(I2),I2=257,384)  
231  C  C
232  C  C
233  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 244  C     5. Sea ice               - assume
244  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
245  C     7. Surface geopotential  - to be done when orography is in  C     7. Surface geopotential  - to be done when orography is in
246  C                                dynamical kernel. Assume 0. for now.  C                                dynamical kernel. Assume 0. for now.
247        mnthIndex = INT(tYear*12.)+1  
248  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  
249         DO J=1,sNy         DO J=1,sNy
250          DO I=1,sNx          DO I=1,sNx
251           I2 = (sNx)*(J-1)+I           I2 = (sNx)*(J-1)+I
252           alb0(I2,myThid) = 0.           alb0(I2,myThid) = 0.
 c        alb0(I2,myThid) = aim_albedo(I,J,bi,bj)/100.  
253           alb0(I2,myThid) = aim_albedo(I,J,bi,bj)           alb0(I2,myThid) = aim_albedo(I,J,bi,bj)
254          ENDDO          ENDDO
255         ENDDO         ENDDO
256  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)  
257         DO J=1,sNy         DO J=1,sNy
258          DO I=1,sNx          DO I=1,sNx
259           I2 = (sNx)*(J-1)+I           I2 = (sNx)*(J-1)+I
# Line 329  C      CLOSE(1) Line 264  C      CLOSE(1)
264          ENDDO          ENDDO
265         ENDDO         ENDDO
266  C  C
267  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 :
268  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.  
269         DO J=1,sNy         DO J=1,sNy
270          DO I=1,sNx          DO I=1,sNx
271           I2 = (sNx)*(J-1)+I           I2 = (sNx)*(J-1)+I
272           soilq1(I2,myThid) = 0.           soilq1(I2,myThid) = 0.
 c        soilq1(I2,myThid) = aim_soilMoisture(I,J,bi,bj)/20.  
273           soilq1(I2,myThid) = aim_soilMoisture(I,J,bi,bj)           soilq1(I2,myThid) = aim_soilMoisture(I,J,bi,bj)
274          ENDDO          ENDDO
275         ENDDO         ENDDO
 C_cnh01      ENDIF  
 C  
 C_cnh01      IF ( FirstCall ) THEN  
276  C      Set snow depth, sea ice to zero for now  C      Set snow depth, sea ice to zero for now
277  C      Land-sea mask ( figure this out from where  C      Land-sea mask ( figure this out from where
278  C                      soil moisture is exactly zero ).  C                      soil moisture is exactly zero ).
# Line 363  C                      soil moisture is Line 288  C                      soil moisture is
288  C      open(77,file='lsmask',form='unformatted')  C      open(77,file='lsmask',form='unformatted')
289  C      write(77) fmask1  C      write(77) fmask1
290  C      close(77)  C      close(77)
291  C_cnh01      ENDIF  
 C  
292  C Addition may 15 . Reset humidity to 0. if negative  C Addition may 15 . Reset humidity to 0. if negative
293  C ---------------------------------------------------  C ---------------------------------------------------
294  #ifdef OLD_AIM_INTERFACE  #ifdef OLD_AIM_INTERFACE
# Line 385  C -------------------------------------- Line 309  C --------------------------------------
309  C     Calculate diagnostics for AIM  C     Calculate diagnostics for AIM
310        CALL AIM_CALC_DIAGS( bi, bj, currentTime, myThid )        CALL AIM_CALC_DIAGS( bi, bj, currentTime, myThid )
311  #endif /* ALLOW_TIMEAVE */  #endif /* ALLOW_TIMEAVE */
 C  
       FirstCall = .FALSE.  
312    
313        CALL AIM_AIM2DYN( bi, bj, currentTime, myThid )        CALL AIM_AIM2DYN( bi, bj, currentTime, myThid )
314  C  
315  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
316    
317        RETURN        RETURN
Legend:
Removed from v.1.8.2.1  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22