/[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.3 by jmc, Tue Mar 6 17:53:02 2001 UTC revision 1.4 by cnh, Tue May 29 19:28:53 2001 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C     $Name$  C     $Name$
3    
4  #include "AIM_OPTIONS.h"  #include "AIM_OPTIONS.h"
5    #undef OLD_AIM_GRIG_MAPPING
6    
7        SUBROUTINE AIM_DO_ATMOS_PHYSICS( phi_hyd, currentTime, myThid )        SUBROUTINE AIM_DO_ATMOS_PHYSICS( phi_hyd,
8         I                                 bi, bj,
9         I                                 currentTime, myThid )
10  C     /==================================================================\  C     /==================================================================\
11  C     | S/R AIM_DO_ATMOS_PHYSICS                                         |  C     | S/R AIM_DO_ATMOS_PHYSICS                                         |
12  C     |==================================================================|  C     |==================================================================|
# Line 15  C     | which can be included as externa Line 18  C     | which can be included as externa
18  C     | tendency routines. Packages should communicate this information  |  C     | tendency routines. Packages should communicate this information  |
19  C     | through common blocks.                                           |  C     | through common blocks.                                           |
20  C     \==================================================================/  C     \==================================================================/
21          IMPLICIT rEAL*8 (A-H,O-Z)
22    
23  C     -------------- Global variables ------------------------------------  C     -------------- Global variables ------------------------------------
24  C     Physics package  C     Physics package
# Line 25  C     Physics package Line 29  C     Physics package
29        INTEGER NLAT        INTEGER NLAT
30        INTEGER NLEV        INTEGER NLEV
31        PARAMETER ( NLON=IX, NLAT=IL, NLEV=KX, NGP=NLON*NLAT )        PARAMETER ( NLON=IX, NLAT=IL, NLEV=KX, NGP=NLON*NLAT )
32  #include "com_physvar.h"  
 #include "com_forcing1.h"  
 #include "Lev_def.h"  
33  C     MITgcm  C     MITgcm
34  #include "EEPARAMS.h"  #include "EEPARAMS.h"
35  #include "PARAMS.h"  #include "PARAMS.h"
36  #include "DYNVARS.h"  #include "DYNVARS.h"
37  #include "GRID.h"  #include "GRID.h"
38    #include "SURFACE.h"
39    #include "AIM_FFIELDS.h"
40    
41    C     Physics package
42    #include "com_physvar.h"
43    #include "com_forcing1.h"
44    #include "Lev_def.h"
45    
46  C     -------------- Routine arguments -----------------------------------  C     -------------- Routine arguments -----------------------------------
47        _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
48        _RL currentTime        _RL currentTime
49          INTEGER myThid
50          INTEGER bi, bj
51    
52  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
53  C     -------------- Local variables -------------------------------------  C     -------------- Local variables -------------------------------------
# Line 103  C     Katm          - Atmospheric K inde Line 114  C     Katm          - Atmospheric K inde
114        real pvoltotNiv5        real pvoltotNiv5
115        SAVE pvoltotNiv5        SAVE pvoltotNiv5
116        real ptotalNiv5        real ptotalNiv5
       INTEGER bi, bj  
117        INTEGER Katm        INTEGER Katm
118    
119  C  C
120        pGround = 1.D5        pGround = 1.D5
121        CPAIR   = 1004        CPAIR   = 1004
122        RD      =  287        RD      =  287
123    
124          CALL AIM_DYN2AIM( bi, bj, currentTime, myThid )
125    
126  C     Assume only one tile per proc. for now  C     Assume only one tile per proc. for now
127        bi  = 1        IG0 = myXGlobalLo+(bi-1)*sNx
128        bj  = 1        JG0 = myYGlobalLo+(bj-1)*sNy
       IG0 = myXGlobalLo  
       JG0 = myYGlobalLo  
129    
130  C        C      
131  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.
# Line 140  C       the mean heave of the base of th Line 151  C       the mean heave of the base of th
151         DO J=1,sNy         DO J=1,sNy
152          DO I=1,sNx          DO I=1,sNx
153           phiTotal(I,J,K) = phiTotal(I,J,K) +           phiTotal(I,J,K) = phiTotal(I,J,K) +
154       &   recip_rhoConst*(phi_hyd(i,j,k,bi,bj))       &   recip_rhoConst*(phi_hyd(i,j,k))
155          ENDDO          ENDDO
156         ENDDO         ENDDO
157        ENDDO        ENDDO
# Line 158  C     _GLOBAL_SUM_R8( phiTSum, myThid ) Line 169  C     _GLOBAL_SUM_R8( phiTSum, myThid )
169  C     ptotalniv5=phiTSum/phiTCount  C     ptotalniv5=phiTSum/phiTCount
170        ptotalniv5=0.        ptotalniv5=0.
171    
172  C     Note the mapping here is only valid for one tile  c_jmc: Because AIM physics LSC is not applied in the stratosphere (top level),
173  C     per proc.  c      ==> move water wapor from the stratos to the surface level.
174          DO J = 1-Oly, sNy+Oly
175           DO I = 1-Olx, sNx+Olx
176    c       k = k_surf(i,j,bi,bj)
177    c       salt(I,J,k,bi,bj) = salt(I,J,k,bi,bj)
178    c    &   + maskC(i,j,Nr,bi,bj)*salt(I,J,Nr,bi,bj)*drF(Nr)*recip_drF(k)
179            salt(I,J,Nr,bi,bj) = 0.
180           ENDDO
181          ENDDO
182    
183    C     Note the mapping here is only valid for one tile per proc.
184        DO K = 1, Nr        DO K = 1, Nr
185         DO J = 1, sNy         DO J = 1, sNy
186          DO I = 1, sNx          DO I = 1, sNx
187           I2 = (sNx)*(J-1)+I           I2 = (sNx)*(J-1)+I
188           Katm = _KD2KA( K )           Katm = _KD2KA( K )
189           UG1(I2,Katm)   = 0.5*(uVel(I,J,K,bi,bj)+uVel(I+1,J,K,bi,bj))           UG1(I2,Katm,myThid)   =
190           VG1(I2,Katm)   = 0.5*(vVel(I,J,K,bi,bj)+vVel(I,J+1,K,bi,bj))       &    0.5*(uVel(I,J,K,bi,bj)+uVel(I+1,J,K,bi,bj))
191  C        Phyiscs works with temperature - not potential temp.           VG1(I2,Katm,myThid)   =
192           TG1(I2,Katm)   = theta(I,J,K,bi,bj)/((pGround/pSurfs(K))**(RD/CPAIR))       &    0.5*(vVel(I,J,K,bi,bj)+vVel(I,J+1,K,bi,bj))
193           QG1(I2,Katm)   = salt(I,J,K,bi,bj)  C        Physics works with temperature - not potential temp.
194           PHIG1(I2,Katm) = (phiTotal(I,J,K)- ptotalniv5 ) + gravity*Hinitial(k)           TG1(I2,Katm,myThid)   = theta(I,J,K,bi,bj)
195         &                  / ((pGround/pSurfs(K))**(RD/CPAIR))
196    c_jmc    QG1(I2,Katm,myThid)   = salt(I,J,K,bi,bj)
197             QG1(I2,Katm,myThid)   = MAX(salt(I,J,K,bi,bj), 0. _d 0)
198             PHIG1(I2,Katm,myThid) = (phiTotal(I,J,K)- ptotalniv5 )
199         &                  + gravity*Hinitial(k)
200    C *NOTE* Fix me for lopped cells
201           if(hFacC(i,j,k,bi,bj).eq.1.) then           if(hFacC(i,j,k,bi,bj).eq.1.) then
202             RHOG1(I2,Katm) = pSurfs(K)/RD/TG1(I2,Katm)             RHOG1(I2,Katm) = pSurfs(K)/RD/TG1(I2,Katm,myThid)
203           else           else
204             RHOG1(I2,Katm)=0.             RHOG1(I2,Katm)=0.
205           endif           endif
206    C *NOTE* Fix me for lopped cells
207          ENDDO          ENDDO
208         ENDDO         ENDDO
209        ENDDO        ENDDO
210    
211    c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
212    c_jmc: add square of surface wind speed (center of C grid) = 2 * KE_surf
213          DO J = 1, sNy
214           DO I = 1, sNx
215             I2 = I+(J-1)*sNx
216             K = k_surf(i,j,bi,bj)
217             Vsurfsq(I2,myThid) = 0.5 * (
218         &                 uVel(I,J,K,bi,bj)*uVel(I,J,K,bi,bj)
219         &               + uVel(I+1,J,K,bi,bj)*uVel(I+1,J,K,bi,bj)
220         &               + vVel(I,J,K,bi,bj)*vVel(I,J,K,bi,bj)
221         &               + vVel(I,J+1,K,bi,bj)*vVel(I,J+1,K,bi,bj)
222         &                        )
223    #ifdef OLD_AIM_GRIG_MAPPING
224    c - to reproduce old results :
225             Katm = _KD2KA( K )
226             Vsurfsq(I2,myThid) =  
227         &     UG1(I2,Katm,myThid)*UG1(I2,Katm,myThid)
228         &   + VG1(I2,Katm,myThid)*VG1(I2,Katm,myThid)
229    #endif /* OLD_AIM_GRIG_MAPPING */
230           ENDDO
231          ENDDO
232    c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
233    
234  C        C      
235  C     Set geopotential surfaces  C     Set geopotential surfaces
236  C     -------------------------  C     -------------------------
237        DO J=1,sNy        DO J=1,sNy
238         DO I=1,sNx         DO I=1,sNx
239          I2 = (sNx)*(J-1)+I          I2 = (sNx)*(J-1)+I
240          IF ( Nlevxy(I2) .NE. 0 ) THEN          IF ( Nlevxy(I2,myThid) .NE. 0 ) THEN
241           PHI0(I2) = gravity*Hinitialw(Nlevxy(I2))           PHI0(I2,myThid) = gravity*Hinitialw(Nlevxy(I2,myThid))
242          ELSE          ELSE
243           PHI0(I2) = 0.           PHI0(I2,myThid) = 0.
244          ENDIF          ENDIF
245         ENDDO         ENDDO
246        ENDDO        ENDDO
247    
248  C  C
249  C     Physics package works with log of surface pressure  C     Physics package works with log of surface pressure
250  C     Get surface pressure from pbot-dpref/dz*Z'  C     Get surface pressure from pbot-dpref/dz*Z'
251        DO J=1,sNy        DO J=1,sNy
252         DO I=1,sNx         DO I=1,sNx
253          I2 = (sNx)*(J-1)+I          I2 = (sNx)*(J-1)+I
254          IF ( Nlevxy(I2) .NE. 0 ) THEN          IF ( Nlevxy(I2,myThid) .NE. 0 ) THEN
255           PNLEVW(I2) = PsurfW(Nlevxy(I2))/pGround           PNLEVW(I2,myThid) = PsurfW(Nlevxy(I2,myThid))/pGround
256          ELSE          ELSE
257  C        Dummy value for land  C        Dummy value for land
258           PNLEVW(I2) = PsurfW(1)/pGround           PNLEVW(I2,myThid) = PsurfW(1)/pGround
259          ENDIF          ENDIF
260          PSLG1(I2) = 0.          PSLG1(I2,myThid) = 0.
261         ENDDO         ENDDO
262        ENDDO        ENDDO
263  cch      write(0,*)  '(PNLEVW(I2),I2=257,384)'  cch      write(0,*)  '(PNLEVW(I2),I2=257,384)'
# Line 215  C Line 267  C
267  C     Physics package needs to know time of year as a fraction  C     Physics package needs to know time of year as a fraction
268        tYear = currentTime/(86400.*360.) -        tYear = currentTime/(86400.*360.) -
269       &        FLOAT(INT(currentTime/(86400.*360.)))       &        FLOAT(INT(currentTime/(86400.*360.)))
270    
271  C  C
272  C     Load external data needed by physics package  C     Load external data needed by physics package
273  C     1. Albedo  C     1. Albedo
# Line 226  C     6. Land sea mask         - infer f Line 279  C     6. Land sea mask         - infer f
279  C     7. Surface geopotential  - to be done when orography is in  C     7. Surface geopotential  - to be done when orography is in
280  C                                dynamical kernel. Assume 0. for now.  C                                dynamical kernel. Assume 0. for now.
281        mnthIndex = INT(tYear*12.)+1        mnthIndex = INT(tYear*12.)+1
282        IF ( mnthIndex .NE. prevMnthIndex .OR.  C_cnh01      IF ( mnthIndex .NE. prevMnthIndex .OR.
283       &     FirstCall ) THEN  C_cnh01     &     FirstCall ) THEN
284         prevMnthIndex = mnthIndex  C_cnh01       prevMnthIndex = mnthIndex
285  C      Read in surface albedo data (input is in % 0-100 )  C      Read in surface albedo data (input is in % 0-100 )
286  C      scale to give fraction between 0-1 for Francos package.  C      scale to give fraction between 0-1 for Francos package.
287  CequChan       WRITE(fNam,'(A,A,A)' ) 'salb.',mnthNam(mnthIndex),'.sun.b'  C      WRITE(fNam,'(A,A,A)' ) 'salb.',mnthNam(mnthIndex),'.sun.b'
288  CequChan       OPEN(1,FILE=fNam(1:14),STATUS='old',FORM='unformatted')  C      OPEN(1,FILE=fNam(1:14),STATUS='old',FORM='unformatted')
289  CequChan       READ(1) tmp4  C      READ(1) tmp4
290  CequChan       CLOSE(1)  C      CLOSE(1)
291  CequChan       DO J=1,nYio         DO J=1,nYio
292  CequChan        DO I=1,nXio          DO I=1,nXio
293  CequChan         tmp4(I,J) = tmp4(I,J)/100.           tmp4(I,J) = aim_albedo(I,J)/100.
294  CequChan        ENDDO          ENDDO
295  CequChan       ENDDO         ENDDO
296         DO J=1,sNy         DO J=1,sNy
297          DO I=1,sNx          DO I=1,sNx
298           I2 = (sNx)*(J-1)+I           I2 = (sNx)*(J-1)+I
299           alb0(I2) = 0.           alb0(I2,myThid) = 0.
300  CequChan         IF ( IG0+I-1 .LE. nxIo .AND. JG0+J-1 .LE. nyIo ) THEN           IF ( IG0+I-1 .LE. nxIo .AND. JG0+J-1 .LE. nyIo ) THEN
301  CequChan          alb0(I2) = tmp4(IG0+I-1,JG0+J-1)            alb0(I2,myThid) = tmp4(IG0+I-1,JG0+J-1)
302  CequChan         ENDIF           ENDIF
303          ENDDO          ENDDO
304         ENDDO         ENDDO
305  C      Read in surface temperature data (input is in absolute temperature)  C      Read in surface temperature data (input is in absolute temperature)
306  CequChan       WRITE(fNam,'(A,A,A)' ) 'tsurf.',mnthNam(mnthIndex),'.sun.b'  C      WRITE(fNam,'(A,A,A)' ) 'tsurf.',mnthNam(mnthIndex),'.sun.b'
307  CequChan       OPEN(1,FILE=fNam(1:15),STATUS='old',FORM='unformatted')  C      OPEN(1,FILE=fNam(1:15),STATUS='old',FORM='unformatted')
308  CequChan       READ(1) tmp4  C      READ(1) tmp4
309  CequChan       CLOSE(1)  C      CLOSE(1)
310         DO J=1,sNy         DO J=1,sNy
311          DO I=1,sNx          DO I=1,sNx
312           I2 = (sNx)*(J-1)+I           I2 = (sNx)*(J-1)+I
313           sst1(I2) = 300.           sst1(I2,myThid) = 300.
314           stl1(I2) = 300.           stl1(I2,myThid) = 300.
315  CequChan         IF ( IG0+I-1 .LE. nxIo .AND. JG0+J-1 .LE. nyIo ) THEN           IF ( IG0+I-1 .LE. nxIo .AND. JG0+J-1 .LE. nyIo ) THEN
316  CequChan          sst1(I2) = tmp4(IG0+I-1,JG0+J-1)            sst1(I2,myThid) = aim_surfTemp(IG0+I-1,JG0+J-1)
317  CequChan          stl1(I2) = tmp4(IG0+I-1,JG0+J-1)            stl1(I2,myThid) = aim_surfTemp(IG0+I-1,JG0+J-1)
318  CequChan         ENDIF           ENDIF
 caja     IF ( I .GE. 64-10 .AND. I .LE. 65+10 ) THEN  
 caja      sst1(I2) = 310.  
 caja      stl1(I2) = 310.  
 caja     ENDIF  
 caja     IF ( I .GE. 64-10 .AND. I .LE. 65+10 ) THEN  
 caja      sst1(I2) = 300.+10.*exp( -((float(I)-64.5)/5.)**2 )  
 caja      stl1(I2) = sst1(I2)  
 caja     ENDIF  
 c_jmc: should not be part of the AIM package :  
          sst1(I2) = 300.+10.*exp( -((float(I)-64.5)/25.)**2 )  
          stl1(I2) = sst1(I2)  
319          ENDDO          ENDDO
320         ENDDO         ENDDO
321  C  C
322  C      Read in soil moisture data (input is in cm in bucket of depth 20cm. )  C      Read in soil moisture data (input is in cm in bucket of depth 20cm. )
323  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.
324  CequChan       WRITE(fNam,'(A,A,A)' ) 'smoist.',mnthNam(mnthIndex),'.sun.b'  C      WRITE(fNam,'(A,A,A)' ) 'smoist.',mnthNam(mnthIndex),'.sun.b'
325  CequChan       OPEN(1,FILE=fNam(1:16),STATUS='old',FORM='unformatted')  C      OPEN(1,FILE=fNam(1:16),STATUS='old',FORM='unformatted')
326  CequChan       READ(1) tmp4  C      READ(1) tmp4
327  CequChan       CLOSE(1)  C      CLOSE(1)
328  CequChan       WRITE(0,*) ' Read file ', fNam(1:16), IG0, JG0  C      WRITE(0,*) ' Read file ', fNam(1:16), IG0, JG0
329  cdj       tmp4 = (tmp4*7.5/20.)*10.  cdj       tmp4 = (tmp4*7.5/20.)*10.
330         DO J=1,sNy         DO J=1,sNy
331          DO I=1,sNx          DO I=1,sNx
332           I2 = (sNx)*(J-1)+I           I2 = (sNx)*(J-1)+I
333           soilq1(I2) = 0.           soilq1(I2,myThid) = 0.
334  CequChan         IF ( IG0+I-1 .LE. nxIo .AND. JG0+J-1 .LE. nyIo ) THEN           IF ( IG0+I-1 .LE. nxIo .AND. JG0+J-1 .LE. nyIo ) THEN
335  CequChan          soilq1(I2) = tmp4(IG0+I-1,JG0+J-1)            soilq1(I2,myThid) = aim_soilMoisture(IG0+I-1,JG0+J-1)/20.
336  CequChan         ENDIF           ENDIF
         ENDDO  
        ENDDO  
 cdj      Soilqmax=MAxval(soilq1)  
        Soilqmax=20.  
 cdj      if(Soilqmax.ne.0.) then  
        DO J=1,sNy  
         DO I=1,sNx  
          I2 = (sNx)*(J-1)+I  
 CequChan         soilq1(I2)=soilq1(I2)/Soilqmax  
          soilq1(I2) = 1.  
337          ENDDO          ENDDO
338         ENDDO         ENDDO
339  cdj      endif  C_cnh01      ENDIF
       ENDIF  
340  C  C
341        IF ( FirstCall ) THEN  C_cnh01      IF ( FirstCall ) THEN
342  C      Set snow depth, sea ice to zero for now  C      Set snow depth, sea ice to zero for now
343  C      Land-sea mask ( figure this out from where soil moisture is exactly zero ).  C      Land-sea mask ( figure this out from where
344    C                      soil moisture is exactly zero ).
345         DO J=1,sNy         DO J=1,sNy
346          DO I=1,sNx          DO I=1,sNx
347           I2 = (sNx)*(J-1)+I           I2 = (sNx)*(J-1)+I
348           fMask1(I2) = 1.           fMask1(I2,myThid) = 1.
349           IF ( soilq1(I2) .EQ. 0. ) fMask1(I2) = 0.           IF ( soilq1(I2,myThid) .EQ. 0. ) fMask1(I2,myThid) = 0.
350           oice1(I2) = 0.           oice1(I2,myThid) = 0.
351           snow1(I2) = 0.           snow1(I2,myThid) = 0.
352          ENDDO          ENDDO
353         ENDDO         ENDDO
354  C      open(77,file='lsmask',form='unformatted')  C      open(77,file='lsmask',form='unformatted')
355  C      write(77) fmask1  C      write(77) fmask1
356  C      close(77)  C      close(77)
357        ENDIF  C_cnh01      ENDIF
358  C  C
359  C Addition may 15 . Reset humidity to 0. if negative  C Addition may 15 . Reset humidity to 0. if negative
360  C ---------------------------------------------------  C ---------------------------------------------------
361        DO K=1,Nr  Caja  DO K=1,Nr
362         DO J=1-OLy,sNy+OLy  Caja   DO J=1-OLy,sNy+OLy
363          DO I=1-Olx,sNx+OLx  Caja    DO I=1-Olx,sNx+OLx
364           IF ( salt(i,j,k,bi,bj) .LT. 0. .OR. K .EQ. Nr ) THEN  Caja     IF ( salt(i,j,k,bi,bj) .LT. 0. .OR. K .EQ. Nr ) THEN
365            salt(i,j,k,bi,bj) = 0.  Caja      salt(i,j,k,bi,bj) = 0.
366           ENDIF  Caja     ENDIF
367          ENDDO  Caja    ENDDO
368         ENDDO  Caja   ENDDO
369        ENDDO  Caja  ENDDO
370  C  
371        CALL PDRIVER( tYear )        CALL PDRIVER( tYear, myThid )
372    
373  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
374  C     Calculate diagnostics for AIM  C     Calculate diagnostics for AIM
# Line 344  C     Calculate diagnostics for AIM Line 376  C     Calculate diagnostics for AIM
376  #endif /* ALLOW_TIMEAVE */  #endif /* ALLOW_TIMEAVE */
377  C  C
378        FirstCall = .FALSE.        FirstCall = .FALSE.
379    
380          CALL AIM_AIM2DYN( bi, bj, currentTime, myThid )
381  C  C
382  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
383    
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22