File: C:\NOAA\NEMS_11731\src\atmos\gfs\dyn\deldifs.f

1           SUBROUTINE DELDIFS(RTE,WE,QME,XE,YE,TEME,
2          &                   RTO,WO,QMO,XO,YO,TEMO,DELTIM,SL,
3          &                   LS_NODE,COEF00,K_LEVEL,
4          &                   hybrid,gen_coord_hybrid,nislfv)
5     !
6           use gfs_dyn_resol_def
7           use gfs_dyn_layout1
8           use gfs_dyn_coordinate_def				! hmhj
9           use gfs_dyn_deldifs_def
10           use gfs_dyn_physcons, rerth => con_rerth
11          &                    ,  rd => con_rd, cp => con_cp
12           IMPLICIT NONE
13     !
14           logical hybrid, gen_coord_hybrid
15           integer nislfv
16           REAL(KIND=KIND_EVOD) RTE(LEN_TRIE_LS,2,LEVS,NTRAC)
17          &,                     WE(LEN_TRIE_LS,2)
18          &,                    QME(LEN_TRIE_LS,2)
19          &,                     XE(LEN_TRIE_LS,2)
20          &,                     YE(LEN_TRIE_LS,2)
21          &,                   TEME(LEN_TRIE_LS,2)
22          &,                     PE(LEN_TRIE_LS,2)
23     !
24           REAL(KIND=KIND_EVOD) RTO(LEN_TRIO_LS,2,LEVS,NTRAC)
25          &,                     WO(LEN_TRIO_LS,2)
26          &,                    QMO(LEN_TRIO_LS,2)
27          &,                     XO(LEN_TRIO_LS,2)
28          &,                     YO(LEN_TRIO_LS,2)
29          &,                   TEMO(LEN_TRIO_LS,2)
30          &,                     PO(LEN_TRIO_LS,2)
31     !
32           REAL(KIND=KIND_EVOD) DELTIM, SL(LEVS)
33     !
34           INTEGER              LS_NODE(LS_DIM,3)
35     !
36     !CMR  LS_NODE(1,1) ... LS_NODE(LS_MAX_NODE,1) : VALUES OF L
37     !CMR  LS_NODE(1,2) ... LS_NODE(LS_MAX_NODE,2) : VALUES OF JBASEV
38     !CMR  LS_NODE(1,3) ... LS_NODE(LS_MAX_NODE,3) : VALUES OF JBASOD
39     !
40           REAL(KIND=KIND_EVOD) COEF00(LEVS,NTRAC)
41     !
42           INTEGER              K_LEVEL
43     !
44           INTEGER              I,IS,IT,JDEL,JDELH,K,KD,KU
45           INTEGER              L,LOCL,N,N0,ND,NP,NPD
46     !
47           INTEGER              INDEV
48           INTEGER              INDOD
49           integer              indev1,indev2
50           integer              indod1,indod2
51           real(kind=kind_evod), parameter :: rkappa = cp / rd
52           REAL(KIND=KIND_EVOD) DN1,REALVAL,RTNP,DF_DK,FACT
53          &,                    SLRD0,FTRD1,RFACT,RFACTRD,RTRD1,FSHK, tem
54     !
55           REAL(KIND=KIND_EVOD), parameter :: CONS0=0.0, CONS1=1.0, CONS2=2.0
56     !
57           INTEGER              INDLSEV,JBASEV
58           INTEGER              INDLSOD,JBASOD
59     !
60           INCLUDE 'function2'
61     !
62     !     print *,' enter deldifs_fd ' 					! hmhj
63     !......................................................................
64     !
65           IF(K_LEVEL.EQ.0) THEN
66     !
67             CALL COUNTPERF(0,15,0.)
68     !!
69             allocate(RTRD(LEVS),RTHK(LEVS),sf(levs))
70             ALLOCATE ( DNE(LEN_TRIE_LS) )
71             ALLOCATE ( DNO(LEN_TRIO_LS) )
72             ALLOCATE ( BKLY(levs) )        					! hmhj
73             ALLOCATE ( CKLY(levs) )        					! hmhj
74             BKLY(:) = 1.0
75             CKLY(:) = 0.0
76             if (gen_coord_hybrid) then					! hmhj
77               DO  k=1,LEVS							! hmhj
78     ! hmhj ak5, bk5, ck5 in gen_coord_hybrid is the same order as model index
79                 BKLY(k)=0.5*(bk5(k)+bk5(k+1))				! hmhj
80                 CKLY(k)=0.5*(ck5(k)+ck5(k+1))*rkappa/thref(k)	        ! hmhj
81     !           if( me.eq.0 )						! hmhj
82     !    &         print*,'sl bkly ckly  in deldif=',k,sl(k),bkly(k),ckly(k)! hmhj
83               enddo								! hmhj
84             else if (hybrid) then						! hmhj
85               DO  k=1,LEVS
86     ! hmhj   sl(k) go bottom to top but bk(k) go top to bottom
87                 BKLY(k)=0.5*(bk5(levs-k+1)+bk5(levs-k+2))/SL(k)
88     !           if( me.eq.0 ) print*,'sl bkly in deldif=',k,sl(k),bkly(k)
89               enddo
90             endif
91     !
92             IF(JCAP.GT.170) THEN
93     !         RECIPROCAL OF TIME SCALE OF DIFFUSION AT REFERENCE WAVENUMBER NP
94               RTNP=(JCAP/170.)**4*1.1/3600
95               NP=JCAP
96               N0=0             ! MAXIMUM WAVENUMBER FOR ZERO DIFFUSION
97               JDEL=8           ! ORDER OF DIFFUSION (EVEN POWER TO RAISE DEL)
98               FSHK=2.2         ! EXTRA HEIGHT-DEPENDENT DIFFUSION FACTOR PER SCALE HEIGHT
99             ELSEIF(JCAP.EQ.170) THEN
100     !         RECIPROCAL OF TIME SCALE OF DIFFUSION AT REFERENCE WAVENUMBER NP
101               RTNP=4*3.E15/(RERTH**4)*FLOAT(80*81)**2
102               NP=JCAP
103               N0=0.55*JCAP     ! MAXIMUM WAVENUMBER FOR ZERO DIFFUSION
104               JDEL=2           ! ORDER OF DIFFUSION (EVEN POWER TO RAISE DEL)
105               FSHK=1.0         ! EXTRA HEIGHT-DEPENDENT DIFFUSION FACTOR PER SCALE HEIGHT
106             ELSEIF(JCAP.EQ.126) THEN					! hmhj
107     !         BELOW HAS BEEN TESTED IN SIHMA-THETA FOR 2 YEAR CFS RUN	! hmhj
108               RTNP=4*3.E15/(RERTH**4)*FLOAT(80*81)**2			! hmhj
109               NP=JCAP							! hmhj
110               N0=0.0           						! hmhj
111               JDEL=4           						! hmhj
112               FSHK=1.0         						! hmhj
113             ELSE
114     !         RECIPROCAL OF TIME SCALE OF DIFFUSION AT REFERENCE WAVENUMBER NP
115               RTNP=1*3.E15/(RERTH**4)*FLOAT(80*81)**2
116               NP=JCAP
117               N0=0.55*JCAP     ! MAXIMUM WAVENUMBER FOR ZERO DIFFUSION
118               JDEL=2           ! ORDER OF DIFFUSION (EVEN POWER TO RAISE DEL)
119               FSHK=1.0         ! EXTRA HEIGHT-DEPENDENT DIFFUSION FACTOR PER SCALE HEIGHT
120             ENDIF
121     !
122             IF (ME.EQ.0) THEN
123               PRINT 6,RTNP,NP,N0,JDEL
124         6     FORMAT(' HORIZONTAL DIFFUSION PARAMETERS'/
125          &  '   EFFECTIVE ',6PF10.3,' MICROHERTZ AT WAVENUMBER ',I4/
126          &  '   MAXIMUM WAVENUMBER FOR ZERO DIFFUSION ',I4/
127          &  '   ORDER OF DIFFUSION ',I2)
128             ENDIF
129     !
130             SLRD0=0.002        ! SIGMA LEVEL AT WHICH TO BEGIN RAYLEIGH DAMPING
131             RTRD1=1./(5*86400) ! RECIPROCAL OF TIME SCALE PER SCALE HEIGHT
132                                !  ABOVE BEGINNING SIGMA LEVEL FOR RAYLEIGH DAMPING
133     !
134             DO K=1,LEVS
135               IF(SL(K).LT.SLRD0) THEN
136                 if (k .gt. levr) then
137                   RTRD(K)=RTRD1*LOG(SLRD0/SL(K)) ** 2
138                 else
139                   RTRD(K)=RTRD1*LOG(SLRD0/SL(K))
140                 endif
141               ELSE
142                 RTRD(K)=0
143               ENDIF
144               RTHK(K)=(SL(K))**LOG(1/FSHK)
145             ENDDO
146     !
147             JDELH=JDEL/2
148             NPD=MAX(NP-N0,0)
149             REALVAL=NPD*(NPD+1)
150             DN1=CONS2*RTNP/REALVAL**JDELH
151     !
152     !......................................................................
153     !
154             DO LOCL=1,LS_MAX_NODE
155                    L=LS_NODE(LOCL,1)
156               JBASEV=LS_NODE(LOCL,2)
157               INDEV=INDLSEV(L,L)
158               DO N=L,JCAP,2
159                 ND=MAX(N-N0,0)
160                 REALVAL=ND*(ND+1)
161                 DNE(INDEV)=DN1*REALVAL**JDELH
162                 INDEV=INDEV+1
163               ENDDO
164             ENDDO
165     !
166     !......................................................................
167     !
168             DO LOCL=1,LS_MAX_NODE
169                    L=LS_NODE(LOCL,1)
170               JBASEV=LS_NODE(LOCL,2)
171               if (mod(L,2).eq.mod(jcap+1,2)) then
172                 DNE(INDLSEV(JCAP+1,L))=CONS0 ! SET THE EVEN (N-L) TERMS OF THE TOP ROW TO ZERO
173               ENDIF
174             ENDDO
175     !
176     !......................................................................
177     !
178             DO LOCL=1,LS_MAX_NODE
179                    L=LS_NODE(LOCL,1)
180               JBASOD=LS_NODE(LOCL,3)
181               INDOD=INDLSOD(L+1,L)
182               DO N=L+1,JCAP,2
183                 ND=MAX(N-N0,0)
184                 REALVAL=ND*(ND+1)
185                 DNO(INDOD)=DN1*REALVAL**JDELH
186                 INDOD=INDOD+1
187               ENDDO
188             ENDDO
189     !
190     !......................................................................
191     !
192             DO LOCL=1,LS_MAX_NODE
193                    L=LS_NODE(LOCL,1)
194               JBASOD=LS_NODE(LOCL,3)
195               if (mod(L,2).ne.mod(jcap+1,2)) then
196                 DNO(INDLSOD(JCAP+1,L))=CONS0 ! SET THE ODD (N-L) TERMS OF THE TOP ROW TO ZERO
197               ENDIF
198             ENDDO
199     !
200     !......................................................................
201     !
202             DO K=1,LEVS
203               KD=MAX(K-1,1)
204               KU=MIN(K+1,LEVS)
205               SF(K)=SL(K)/(SL(KU)-SL(KD))/SQRT(CONS2)     !CONSTANT
206             ENDDO
207     !
208             CALL COUNTPERF(1,15,0.)
209     !!
210             RETURN
211           ENDIF
212     !
213     !......................................................................
214     !
215           CALL COUNTPERF(0,13,0.)
216     !!
217           K=K_LEVEL
218     !!
219     !
220     !     TEM = COEF00(K,1) * BKLY(K)
221     !
222           DO LOCL=1,LS_MAX_NODE
223                   L=LS_NODE(LOCL,1)
224              JBASEV=LS_NODE(LOCL,2)
225              IF (L.EQ.0) THEN
226                 N0=2
227              ELSE
228                 N0=L
229              ENDIF
230              indev1 = indlsev(N0,L)
231              if (mod(L,2).eq.mod(jcap+1,2)) then
232                indev2 = indlsev(jcap+1,L)
233              else
234                indev2 = indlsev(jcap  ,L)
235              endif
236     !!       DO N = N0, JCAP+1, 2
237              DO INDEV = indev1 , indev2
238      
239                FACT             = DELTIM*DNE(INDEV)*RTHK(K)
240                RFACT            = CONS1/(CONS1+FACT)
241                RFACTRD          = CONS1/(CONS1+FACT+DELTIM*RTRD(K))
242      
243                WE(INDEV,1)      = WE(INDEV,1)*RFACTRD
244                WE(INDEV,2)      = WE(INDEV,2)*RFACTRD
245      
246                XE(INDEV,1)      = XE(INDEV,1)*RFACTRD
247                XE(INDEV,2)      = XE(INDEV,2)*RFACTRD
248      
249                if( nislfv.le.1 ) then	! hmhj
250                RTE(INDEV,1,1,1) = RTE(INDEV,1,1,1)*RFACT
251                RTE(INDEV,2,1,1) = RTE(INDEV,2,1,1)*RFACT
252                endif			! hmhj
253     
254                PE(INDEV,1)=BKLY(K)*QME(INDEV,1)+CKLY(K)*TEME(INDEV,1)	! hmhj
255                PE(INDEV,2)=BKLY(K)*QME(INDEV,2)+CKLY(K)*TEME(INDEV,2)	! hmhj
256      
257                YE(INDEV,1)      =  ( YE(INDEV,1)+
258          X                 FACT*COEF00(K,1)* PE(INDEV,1) )*RFACT		! hmhj
259      
260                YE(INDEV,2)      = ( YE(INDEV,2) +
261          X                 FACT*COEF00(K,1)* PE(INDEV,2) )*RFACT		! hmhj
262      
263              ENDDO
264            ENDDO
265     !
266     !......................................................................
267     !
268     !      DO L = 0, JCAP
269            DO LOCL=1,LS_MAX_NODE
270                   L=LS_NODE(LOCL,1)
271              JBASOD=LS_NODE(LOCL,3)
272              indod1 = indlsod(L+1,L)
273              if (mod(L,2).eq.mod(jcap+1,2)) then
274                indod2 = indlsod(jcap  ,L)
275              else
276                indod2 = indlsod(jcap+1,L)
277              endif
278     !        DO N = L+1, JCAP+1, 2
279              DO INDOD = indod1 , indod2
280      
281                FACT             = DELTIM*DNO(INDOD)*RTHK(K)
282                RFACT            = CONS1/(CONS1+FACT)
283                RFACTRD          = CONS1/(CONS1+FACT+DELTIM*RTRD(K))
284      
285                WO(INDOD,1)      = WO(INDOD,1)*RFACTRD
286                WO(INDOD,2)      = WO(INDOD,2)*RFACTRD
287      
288                XO(INDOD,1)      = XO(INDOD,1)*RFACTRD
289                XO(INDOD,2)      = XO(INDOD,2)*RFACTRD
290     
291                if( nislfv.le.1 ) then	! hmhj
292                RTO(INDOD,1,1,1) = RTO(INDOD,1,1,1)*RFACT
293                RTO(INDOD,2,1,1) = RTO(INDOD,2,1,1)*RFACT
294                endif			! hmhj
295     
296                PO(INDOD,1)=BKLY(K)*QMO(INDOD,1)+CKLY(K)*TEMO(INDOD,1)	! hmhj
297                PO(INDOD,2)=BKLY(K)*QMO(INDOD,2)+CKLY(K)*TEMO(INDOD,2)	! hmhj
298      
299                YO(INDOD,1)      = ( YO(INDOD,1)+
300          X                 FACT*COEF00(K,1)* PO(INDOD,1) )*RFACT		! hmhj
301      
302                YO(INDOD,2)      = ( YO(INDOD,2)+
303          X                 FACT*COEF00(K,1)* PO(INDOD,2) )*RFACT		! hmhj
304      
305              ENDDO
306            ENDDO
307     !
308     !......................................................................
309     !
310     !
311           if( nislfv.le.1 ) then
312     !
313           if (ntrac .gt. 1) then
314             do it=2,ntrac
315     !
316     !          DO L = 0, JCAP
317                DO LOCL=1,LS_MAX_NODE
318                       L=LS_NODE(LOCL,1)
319                  JBASEV=LS_NODE(LOCL,2)
320                  IF (L.EQ.0) THEN
321                     N0=2
322                  ELSE
323                     N0=L
324                  ENDIF
325                  indev1 = indlsev(N0,L)
326                  if (mod(L,2).eq.mod(jcap+1,2)) then
327                    indev2 = indlsev(jcap+1,L)
328                  else
329                    indev2 = indlsev(jcap  ,L)
330                  endif
331     !            DO N = N0, JCAP+1, 2
332                  DO INDEV = indev1 , indev2
333      
334                    FACT              = DELTIM*DNE(INDEV)*RTHK(K)
335                    RFACT             = CONS1/(CONS1+FACT)
336      
337     !              RTE(INDEV,1,1,IT) =   RTE(INDEV,1,1,IT) * RFACT
338     !              RTE(INDEV,2,1,IT) =   RTE(INDEV,2,1,IT) * RFACT
339     
340                    PE(INDEV,1)=BKLY(K)*QME(INDEV,1)+CKLY(K)*TEME(INDEV,1)	! hmhj
341                    PE(INDEV,2)=BKLY(K)*QME(INDEV,2)+CKLY(K)*TEME(INDEV,2)	! hmhj
342     
343                    RTE(INDEV,1,1,IT) = ( RTE(INDEV,1,1,IT)+
344          X                 FACT*COEF00(K,it)* PE(INDEV,1) )*RFACT		! hmhj
345      
346                    RTE(INDEV,2,1,IT) = ( RTE(INDEV,2,1,IT)+
347          X                 FACT*COEF00(K,it)* PE(INDEV,2) )*RFACT		! hmhj
348      
349                  ENDDO
350                ENDDO
351     !
352     !......................................................................
353     !
354     !         DO L = 0, JCAP
355                DO LOCL=1,LS_MAX_NODE
356                       L=LS_NODE(LOCL,1)
357                  JBASOD=LS_NODE(LOCL,3)
358                  indod1 = indlsod(L+1,L)
359                  if (mod(L,2).eq.mod(jcap+1,2)) then
360                    indod2 = indlsod(jcap  ,L)
361                  else
362                    indod2 = indlsod(jcap+1,L)
363                  endif
364     !            DO N = L+1, JCAP+1, 2
365                  DO INDOD = indod1 , indod2
366      
367                    FACT              = DELTIM*DNO(INDOD)*RTHK(K)
368                    RFACT             = CONS1/(CONS1+FACT)
369      
370     !              RTO(INDOD,1,1,IT) =   RTO(INDOD,1,1,IT) * RFACT
371     !              RTO(INDOD,2,1,IT) =   RTO(INDOD,2,1,IT) * RFACT
372     !
373                    PO(INDOD,1)=BKLY(K)*QMO(INDOD,1)+CKLY(K)*TEMO(INDOD,1)	! hmhj
374                    PO(INDOD,2)=BKLY(K)*QMO(INDOD,2)+CKLY(K)*TEMO(INDOD,2)	! hmhj
375     
376                    RTO(INDOD,1,1,IT) = ( RTO(INDOD,1,1,IT)+
377          X                 FACT*COEF00(K,it)* PO(INDOD,1) )*RFACT		! hmhj
378      
379                    RTO(INDOD,2,1,IT) = ( RTO(INDOD,2,1,IT)+
380          X                 FACT*COEF00(K,it)* PO(INDOD,2) )*RFACT		! hmhj
381      
382                  ENDDO
383                ENDDO
384             enddo                       ! Tracer do loop end
385           endif
386     !
387           endif	! nislfv
388     !
389           CALL COUNTPERF(1,13,0.)
390     
391     !     print *,' leave deldifs_fd ' 					! hmhj
392     !!
393           RETURN
394           END
395