/[MITgcm]/MITgcm/pkg/fizhi/fizhi_rayleigh.F
ViewVC logotype

Diff of /MITgcm/pkg/fizhi/fizhi_rayleigh.F

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

revision 1.1 by molod, Wed Oct 20 18:27:36 2004 UTC revision 1.2 by molod, Fri Oct 22 14:52:14 2004 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3                                                                                                                                                                                                
4  #include "FIZHI_OPTIONS.h"  #include "FIZHI_OPTIONS.h"
5        subroutine rayleigh(myid,pres,preskappa,psurf,u,v,im,jm,lm,bi,bj,        subroutine rayleigh(myid,pres,pkap,pekap,zsurf,u,v,t,s,im,jm,lm,
6       .                                                      rfu,rfv,rft)       .                                                bi,bj,rfu,rfv,rft)
7  C **********************************************************************  C **********************************************************************
8  C  C
9  C  PURPOSE  C  PURPOSE
10  C     To implement Rayleigh Friction  C     To implement Rayleigh Friction
11  C  C
12  C  ARGUMENTS   DESCRIPTION  C  ARGUMENTS   DESCRIPTION
13    C
14    C  INPUT:
15  C     MYID .... PROCESS(OR) NUMBER  C     MYID .... PROCESS(OR) NUMBER
16  C     PRES .... VALUE OF MID-LEVEL PRESSURE  C     PRES .... MID-LEVEL PRESSURE IN MB
17  C     PZ ...... VALUE OF SURFACE PRESSURE - PTOP  C     PKAP .... MID-LEVEL PRESSURE ** KAPPA
18  C     UZ ...... VALUE OF U-WIND IN MODEL FOR N-1 TIME STEP  C     PEKAP ... EDGE-LEVEL PRESSURE ** KAPPA
19  C     VZ ...... VALUE OF V-WIND IN MODEL FOR N-1 TIME STEP  C     ZSURF ... SURFACE ELEVATION IN M
20  C     IM ...... NUMBER OF LONGITUDE POINTS IN MODEL  C     U ....... U-WIND
21  C     JM ...... NUMBER OF LATITUDE  POINTS IN MODEL  C     V ....... V-WIND
22  C     LM ...... NUMBER OF VERTICAL  LEVELS IN MODEL  C     TH ...... THETA (ACTUALLY REAL THETA * P0**KAPPA) IN K
23    C     S  ...... SPECIFIC HUMIDITY (KG/KG)
24    C     IM ...... NUMBER OF LONGITUDE POINTS
25    C     JM ...... NUMBER OF LATITUDE  POINTS
26    C     LM ...... NUMBER OF VERTICAL  LEVELS
27  C     BI ...... X-DIRECTION PROCESSOR INDEX  C     BI ...... X-DIRECTION PROCESSOR INDEX
28  C     BJ ...... Y-DIRECTION PROCESSOR INDEX  C     BJ ...... Y-DIRECTION PROCESSOR INDEX
29  C     RFU ..... VALUE OF U-WIND TENDENCY  C  OUTPUT:
30  C     RFV ..... VALUE OF V-WIND TENDENCY  C     RFU ..... U-WIND TENDENCY
31  C     RFT ..... VALUE OF THETA  TENDENCY  C     RFV ..... V-WIND TENDENCY
32    C     RFT ..... THETA  TENDENCY
33  C  C
34  C **********************************************************************  C **********************************************************************
35    
# Line 35  C ************************************** Line 42  C **************************************
42  #endif  #endif
43    
44        integer myid,im,jm,lm,bi,bj        integer myid,im,jm,lm,bi,bj
45          _RL zsurf(im,jm),pres(im,jm,lm),pkap(im,jm,lm),pekap(im,jm,lm+1)
46        real psurf(im,jm)        _RL u(im,jm,lm),v(im,jm,lm),t(im,jm,lm),s(im,jm,lm)
47        real pres(im,jm,lm)        _RL rfu(im,jm,lm),rfv(im,jm,lm),rft(im,jm,lm)
       real preskappa(im,jm,lm)  
       real u(im,jm,lm)  
       real v(im,jm,lm)  
       real rfu(im,jm,lm)  
       real rfv(im,jm,lm)  
       real rft(im,jm,lm)  
48    
49        integer  i,j,L        integer  i,j,L
50        real rf(im,jm,lm)        _RL rf(im,jm,lm)
51        real z(im,jm,lm)        _RL z(im,jm,lm)
52        real z1(im,jm,lm)        _RL dz(im,jm,lm)
53        real z2(im,jm,lm)        _RL cpog, cpinv, virtcon, getcon
       real cpinv, getcon  
54    
55  C **********************************************************************  C **********************************************************************
56  C ****         APPLY RAYLEIGH FRICTION TO WIND TENDENCIES            ***  C ****   APPLY RAYLEIGH FRICTION TO WIND (INCLUDE HEATING)           ***
57  C **********************************************************************  C **********************************************************************
58    
59          cpog = getcon('CP')/getcon('GRAVITY')
60        cpinv = 1.0/getcon('CP')        cpinv = 1.0/getcon('CP')
61          virtcon = getcon('VIRTCON')
62    
63        do L=1,lm        do L=1,lm
64           do j=1,jm
65        do j=1,jm         do i=1,im
66        do i=1,im          dz(i,j,L) = cpog * (pekap(i,j,L+1)-pekap(i,j,L)) * t(i,j,L) *
67         .                (1.+virtcon*s(i,j,L))
68         z1(i,j,L) = -7e2*log(pres(i,j,1)/psurf(i,j))         enddo
69         z2(i,j,L) = -7e3*log(pres(i,j,2)/psurf(i,j))         enddo
        z(i,j,L) = -7e3*log(pres(i,j,L)/psurf(i,j))  
        rf(i,j,L) = 0.40*(1+tanh((z(i,j,L)-z2(i,j,L))/z1(i,j,L)))/86400  
   
        rfu(i,j,L) = - rf(i,j,L) * u(i,j,L)  
        rfv(i,j,L) = - rf(i,j,L) * v(i,j,L)  
        rft(i,j,L) = -(u(i,j,L)*rfu(i,j,L) + v(i,j,L)*rfv(i,j,L) )*cpinv  
      .                        /preskappa(i,j,L)  
   
       enddo  
70        enddo        enddo
71    
72        if( irfu.ne.0 ) then        do L=1,lm
73        do j=1,jm         do j=1,jm
74        do i=1,im         do i=1,im
75        qdiag(i,j,irfu+L-1,bi,bj) = qdiag(i,j,irfu+L-1,bi,bj) +          dz(i,j,L) = cpog * (pekap(i,j,L+1)-pekap(i,j,L)) * t(i,j,L) *
76       .                 rfu(i,j,L)*86400       .                (1.+virtcon*s(i,j,L))
77        enddo         enddo
78           enddo
79        enddo        enddo
       endif  
80    
       if( irfv.ne.0 ) then  
81        do j=1,jm        do j=1,jm
82        do i=1,im        do i=1,im
83        qdiag(i,j,irfv+L-1,bi,bj) = qdiag(i,j,irfv+L-1,bi,bj) +         z(i,j,lm) = zsurf(i,j) +  0.5 * dz(i,j,lm)
      .                 rfv(i,j,L)*86400  
84        enddo        enddo
85        enddo        enddo
       endif  
86    
87        if( irft.ne.0 ) then        do L=lm-1,1,-1
88        do j=1,jm         do j=1,jm
89        do i=1,im         do i=1,im
90        qdiag(i,j,irft+L-1,bi,bj) = qdiag(i,j,irft+L-1,bi,bj) +          z(i,j,L) = z(i,j,L+1) + 0.5 * (dz(i,j,L)+dz(i,j,L+1))
91       .                 rft(i,j,L)*86400         enddo
92           enddo
93        enddo        enddo
94    
95          do L=1,lm
96           do j=1,jm
97           do i=1,im
98            rf(i,j,L) = (2./3.)*(1+tanh((z(i,j,L)-80000.)/5000.))/86400.
99            rfu(i,j,L) = - rf(i,j,L) * u(i,j,L)
100            rfv(i,j,L) = - rf(i,j,L) * v(i,j,L)
101            rft(i,j,L) = -(u(i,j,L)*rfu(i,j,L) + v(i,j,L)*rfv(i,j,L) )*cpinv
102         .                        /pkap(i,j,L)
103           enddo
104           enddo
105        enddo        enddo
       endif  
106    
107    #ifdef ALLOW_DIAGNOSTICS
108          do L=1,lm
109           if( irfu.ne.0 ) then
110           do j=1,jm
111           do i=1,im
112            qdiag(i,j,irfu+L-1,bi,bj) = qdiag(i,j,irfu+L-1,bi,bj) +
113         .                  rfu(i,j,L)*86400
114           enddo
115           enddo
116           endif
117    
118           if( irfv.ne.0 ) then
119           do j=1,jm
120           do i=1,im
121            qdiag(i,j,irfv+L-1,bi,bj) = qdiag(i,j,irfv+L-1,bi,bj) +
122         .                  rfv(i,j,L)*86400
123           enddo
124           enddo
125           endif
126    
127           if( irft.ne.0 ) then
128           do j=1,jm
129           do i=1,im
130            qdiag(i,j,irft+L-1,bi,bj) = qdiag(i,j,irft+L-1,bi,bj) +
131         .                  rft(i,j,L)*86400
132           enddo
133           enddo
134           endif
135        enddo        enddo
136    #endif
137    
138        nrfu = nrfu + 1  #ifdef ALLOW_DIAGNOSTICS
139        nrfv = nrfv + 1        if( (bi.eq.1) .and. (bj.eq.1) ) then
140        nrft = nrft + 1         nrfu = nrfu + 1
141           nrfv = nrfv + 1
142           nrft = nrft + 1
143          endif
144    #endif
145    
146        return        return
147        end        end
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22