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	#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