igsm/src_chem/chemdiff.F

#include "ctrparam.h"

!       ============================================================
!
!      CHEMDIFF.F:    Subroutine for calculating horizontal
!                       diffusion of MIT Global Chemistry Model
!
!       ------------------------------------------------------------
!
!       Author:         Chien Wang
!                       MIT Joint Program on Science and Policy
!                               of Global Change
!
!       ----------------------------------------------------------
! 
!       Revision History:
!       
!       When    Who             What
!       ----    ----------      -------
!       013096  Chien Wang      rev.    
!       080100  Chien Wang      repack based on CliChem3 & add cpp
!       051804  Chien Wang      rev.
!
!       ==========================================================

      Subroutine chemdiff(ifdiff,x00,x11,dta)

#include "chem_para"
#include "chem_com"
#include "BD2G04.COM"

      dimension x00  (nlon,nlat,nlev)
      dimension x11  (nlon,nlat,nlev)

      dimension dcdy(nlat,nlev)

#if ( defined CPL_CHEM )

c-------------------------------------------------------
c   Definitions of parameters:
c
      istart=1
      iend  =nlon
        
c
c 013096
c   fktdif span from 2.e6 in the first three years to
c       5.e5 or 1.e6 in twenty years and maintain this value
c       thereafter:
c
c       xxx    = float(myyear - 3)/20.0
c       xxx    = amin1(1.0,amax1(0.0,xxx))
c       fktdif = (20.0 - xxx * 10.0)*1.e5       ! m2/s

c      fktdif = 4.e6     !m2/s

c       if(ifdiff.eq.1)then
c         fktdif = 2.e6
c       else if(ifdiff.eq.2)then
c         fktdif = 3.e6
c       endif

c 111596:
c       fktdif = float(ifdiff)*1.e6
        fktdif = float(ifdiff)*1.e5

c=====
c   Calculate dcdy:
c
      do i=istart,iend
      do j=2,nlat
      do k=1,nlev
        dcdy(j,k)=(x11(i,j,k)-x11(i,j-1,k))
     &           /dyv(j)
      end do
      end do
      end do

c=====
c   Calculate meridional eddy diffusion:
c
      do k=1,nlev
        paver  = 0.5*(p00(1,1)+p00(1,2))
        fluxl  =-fktdif
     &         /dyv(2)*dcdy(2,k)*dta 
     &         * paver
        fluxl=max(-0.5*x00(1,2,k), min(0.5*x00(1,1,k),fluxl))
      do j=2,nlat1
        paver  = 0.5*(p00(1,j)+p00(1,j+1))
        fluxr     =-fktdif
     &            /dyv(j+1)*dcdy(j+1,k)*dta
     &         * paver
        fluxr=max(-0.5*x00(1,j+1,k),min(0.5*x00(1,j,k),fluxr))
        x00(1,j,k)=x00(1,j,k)-(fluxr-fluxl)
        fluxl=fluxr
       end do
       end do

c       call chemcheck(x00)

#endif

        return
         end

1	jscott	1.1	#include "ctrparam.h"
2
3			! ============================================================
4			!
5			! CHEMDIFF.F: Subroutine for calculating horizontal
6			! diffusion of MIT Global Chemistry Model
7			!
8			! ------------------------------------------------------------
9			!
10			! Author: Chien Wang
11			! MIT Joint Program on Science and Policy
12			! of Global Change
13			!
14			! ----------------------------------------------------------
15			!
16			! Revision History:
17			!
18			! When Who What
19			! ---- ---------- -------
20			! 013096 Chien Wang rev.
21			! 080100 Chien Wang repack based on CliChem3 & add cpp
22			! 051804 Chien Wang rev.
23			!
24			! ==========================================================
25
26			Subroutine chemdiff(ifdiff,x00,x11,dta)
27
28			#include "chem_para"
29			#include "chem_com"
30			#include "BD2G04.COM"
31
32			dimension x00 (nlon,nlat,nlev)
33			dimension x11 (nlon,nlat,nlev)
34
35			dimension dcdy(nlat,nlev)
36
37			#if ( defined CPL_CHEM )
38
39			c-------------------------------------------------------
40			c Definitions of parameters:
41			c
42			istart=1
43			iend =nlon
44
45			c
46			c 013096
47			c fktdif span from 2.e6 in the first three years to
48			c 5.e5 or 1.e6 in twenty years and maintain this value
49			c thereafter:
50			c
51			c xxx = float(myyear - 3)/20.0
52			c xxx = amin1(1.0,amax1(0.0,xxx))
53			c fktdif = (20.0 - xxx * 10.0)*1.e5 ! m2/s
54
55			c fktdif = 4.e6 !m2/s
56
57			c if(ifdiff.eq.1)then
58			c fktdif = 2.e6
59			c else if(ifdiff.eq.2)then
60			c fktdif = 3.e6
61			c endif
62
63			c 111596:
64			c fktdif = float(ifdiff)*1.e6
65			fktdif = float(ifdiff)*1.e5
66
67			c=====
68			c Calculate dcdy:
69			c
70			do i=istart,iend
71			do j=2,nlat
72			do k=1,nlev
73			dcdy(j,k)=(x11(i,j,k)-x11(i,j-1,k))
74			& /dyv(j)
75			end do
76			end do
77			end do
78
79			c=====
80			c Calculate meridional eddy diffusion:
81			c
82			do k=1,nlev
83			paver = 0.5*(p00(1,1)+p00(1,2))
84			fluxl =-fktdif
85			& /dyv(2)dcdy(2,k)dta
86			& * paver
87			fluxl=max(-0.5x00(1,2,k), min(0.5x00(1,1,k),fluxl))
88			do j=2,nlat1
89			paver = 0.5*(p00(1,j)+p00(1,j+1))
90			fluxr =-fktdif
91			& /dyv(j+1)dcdy(j+1,k)dta
92			& * paver
93			fluxr=max(-0.5x00(1,j+1,k),min(0.5x00(1,j,k),fluxr))
94			x00(1,j,k)=x00(1,j,k)-(fluxr-fluxl)
95			fluxl=fluxr
96			end do
97			end do
98
99			c call chemcheck(x00)
100
101			#endif
102
103			return
104			end
105