/[MITgcm]/MITgcm_contrib/jscott/igsm/src/forcedozone.F
ViewVC logotype

Contents of /MITgcm_contrib/jscott/igsm/src/forcedozone.F

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


Revision 1.1 - (show annotations) (download)
Fri Aug 11 19:35:30 2006 UTC (18 years, 11 months ago) by jscott
Branch: MAIN
CVS Tags: HEAD
atm2d package

1
2 #include "ctrparam.h"
3
4 ! ==========================================================
5 !
6 ! FORCEDOZONE.F: Routine to get ozone deviations from
7 ! specified file and interpolate
8 ! to the appropriate data
9 !
10 ! ----------------------------------------------------------
11 !
12 ! Revision History:
13 !
14 ! When Who What
15 ! ---- ---------- -------
16 ! 073100 Chien Wang repack based on CliChem3 & M24x11,
17 ! and add cpp.
18 !
19 ! 110100 Added data for Jan. 2000,
20 ! dimension of o3f was changed from nm to nm+1
21 !
22 ! 011503 Chris Forest Converted to use GISS 1850-2050 ozone data
23 ! uses nlev = 14 to adjust ozone profile in three
24 ! layers above 10 hPa top of dynamics levels
25 ! (layer edges at 5, 2, and 0.00001 hPa)
26 !
27 ! ==========================================================
28
29 subroutine getforcedozone( jyear, jday )
30
31 parameter (nm=1728)
32 c parameter (nm =252)
33
34 #include "chem_para"
35 #include "chem_com"
36
37 common /o3data/o3_data
38 character * 120 o3_data
39 C
40 c Routine to get ozone deviations from specified file and interpolate
41 c to the appropriate date.
42 C
43 C Author: Chris E Forest.
44 C Added Mar 18, 1998
45 C
46 C tropmass = 28.97296245*1.e-3*0.8/P0
47 C pxxx = dpl(l)*part
48 C
49 C ULGAS(L,3)=o3(ILON,jyyy,l)/48.0
50 C & *pxxx*tropmass
51 C
52 C o3 = ppb(m)
53 C
54 C 48 = mol weight of o3
55 C
56 C ULGAS = cm^3 (STP)/cm^2
57
58 logical first
59 integer year0, nlevp3 ,yrmax
60 ! real date, o3f(nlat,nlev+3,nm), o3dummy ! 03/20/2006
61 real date, o3f(nlat,nlev+3,nm+12), o3dummy
62 c real date, o3f(24,11,nm+1), o3dummy
63 character*43 filen
64 DIMENSION JDOFM(13)
65 data first / .true. /
66 DATA JDOFM /0,31,59,90,120,151,181,212,243,273,304,334,365/
67 data year0 / 1859 /
68
69 c year0 = 1974 for o3hadley.dat
70 c year0 = 1979 for o3giss.MR.dat
71 c year0 = 1850 for o3.1850_2050.giss.46x11.dat
72 c year0 = 1859 for o3.1859_2002.giss.46x11.dat
73 c
74 C filen = '/data11/ceforest/forcing/dummyozone.dat'
75 c filen = '/data11/ceforest/forcing/o3giss.MR.dat'
76 C filen = '/data11/ceforest/forcing/o3hadley.dat'
77
78 if (first) then
79 print *,'GetForcedOzone: reading data for 14 layers'
80 print *,'nlat=',nlat,' nlev=',nlev
81 open(unit=500, file=o3_data,status='old')
82 c
83 c GISS: nm = 252 , HadCM2: nm=265, 1850-2050: nm=2412
84 c
85 c nlev = nlev + 3 to include 3 layers above 10 hPa.
86 c
87
88 nlevp3 = nlev + 3
89
90 do i=1,nm
91 do k=1, nlevp3
92 do j = 1,nlat
93 read(500,'(E12.4)') o3dummy
94 o3f(j,k,i) = o3dummy*1.e9 ! convert to ppb(m)
95 enddo
96 enddo
97 enddo
98 close(500)
99 C Added data for Jan. 2000 to allow simulation until Jan 1 2000
100 c NOT necessary with 1850-2050 GISS ozone data
101 c do k=1, nlev
102 c do j = 1,nlat
103 c o3f(j,k,nm+1) = o3f(j,k,nm)
104 c enddo
105 c enddo
106 ! 03/20/2006 Data for anaditional year are added
107 ! really only data for January are used
108 do m=1, 12
109 do k=1, nlev
110 do j = 1,nlat
111 o3f(j,k,nm+m) = o3f(j,k,nm-12+m)
112 enddo
113 enddo
114 enddo
115 first=.false.
116 nyears=nm/12
117 yrmax=year0+nyears-1
118 print *,'Ozone data for ',nyears,' years '
119 print *,'from ',year0,' to ',yrmax
120 endif
121
122
123 ! iy = jyear - year0 ! 03/17/06
124 jyear1=min(jyear,yrmax)
125 iy = jyear1 - year0
126 ! print *,'Ozone for year ',jyear1,iy
127 if (iy.lt.0) then
128 do i=1,nlat
129 do j=1,nlevp3
130 o3dev(1,i,j) = 0.
131 enddo
132 enddo
133 else
134 C--- o3interp(o3f, o3dev, date, year0)
135
136 jdd = jday-15
137
138 if (iy .eq. 0 .and. jdd.le.0) then
139 do i=1,nlat
140 do j=1,nlevp3
141 dl = ( 31. + float(jdd) )/31.
142 o3dev(1,i,j) = o3f(i,j,1)*dl
143 enddo
144 enddo
145 else
146 im =12
147 do i=1,12
148 if (jdd.gt.jdofm(i) .and. jdd.le.jdofm(i+1)) then
149 im = i
150 endif
151 enddo
152 if (iy .ge. 1 .and. im.eq.12 .and. jdd.le.0 ) iy = iy-1
153 C if (jdd.le.0) iy = iy-1
154 C endif
155
156 dl = ( float(jdd) - jdofm(im))/(jdofm(im+1) - jdofm(im))
157 if (im.eq.12.and.jdd.le.0) dl = ( 31. + float(jdd) )/31.
158
159 imm = im + 12*iy
160 ! print *,jyear1,jday,iy,im,imm
161
162 if(imm.lt.1.or.imm.ge.nm+1)then
163 print *,' error in o3interp im=',imm
164 stop 25
165 endif
166
167 do j=1,nlat
168 do k=1,nlevp3
169 o3dev(1,j,k) = o3f(j,k,imm)*(1.-dl)+o3f(j,k,imm+1)*dl
170 enddo
171 enddo
172 C print *,jyear,jday,iy,im,dl,o3dev(0,0,9)
173 endif
174 endif
175
176 c print *,'GetForcedOzone: JYEAR, JDAY, IY=', jyear,jday, iy
177 c if (iy.ge.0) print *,'GetForcedOzone: imm=', imm
178
179 return
180 end
181
182

  ViewVC Help
Powered by ViewVC 1.1.22