/[MITgcm]/MITgcm/verification/fizhi-gridalt-hs/code/do_fizhi.F
ViewVC logotype

Contents of /MITgcm/verification/fizhi-gridalt-hs/code/do_fizhi.F

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


Revision 1.6 - (show annotations) (download)
Thu Oct 21 15:19:18 2004 UTC (19 years, 5 months ago) by molod
Branch: MAIN
Changes since 1.5: +3 -2 lines
Accomodate new do_fizhi argument list

1 C $Header: /u/gcmpack/MITgcm/verification/fizhi-gridalt-hs/code/do_fizhi.F,v 1.5 2004/10/15 21:15:06 molod Exp $
2 C $Name: $
3 #include "FIZHI_OPTIONS.h"
4 subroutine do_fizhi(myid,
5 . idim1,idim2,jdim1,jdim2,Nrphin,Nsxin,Nsyin,im1,im2,jm1,jm2,bi,bj,
6 . nchp,nchptot,nchpland,
7 . uphy,vphy,thphy,sphy,pephy,lons,lats,phi0surf,
8 . ctmt,xxmt,yymt,zetamt,xlmt,khmt,tke,
9 . tgz,sst,sice,phis_var,landtype,fracland,emiss,albnirdr,albnirdf,
10 . albvisdr,albvisdf,ityp,chfr,alai,agrn,igrd,chlat,chlon,
11 . tcanopy,tdeep,ecanopy,swetshal,swetroot,swetdeep,snodep,capac,
12 . o3,qstr,co2,cfc11,cfc12,cfc22,n2o,methane,
13 . iras,nlwcld,cldtotlwin,cldraslwin,cldlsplwin,nlwlz,lwlzin,
14 . nswcld,cldtotswin,cldrasswin,cldlspswin,nswlz,swlzin,imstturbsw,
15 . imstturblw,qliqaveswin,qliqavelwin,fccaveswin,fccavelwin,
16 . rainconin,rainlspin,snowfallin,
17 . duphy,dvphy,dthphy,dsphy)
18 c-----------------------------------------------------------------------
19 c Interface routine to calculate physics increments - calls fizhi_driver.
20 c Purpose of this routine is to set up arrays local to fizhi and 'save'
21 c them from one iteration to the next, and act as interface between the
22 c model common blocks (held in fizhi_wrapper) and fizhi_driver.
23 c Copies of variables that are 'shadowed' are made here without shadows
24 c for passing to fizhi_driver.
25 c Note: routine is called from inside a bi-bj loop
26 c
27 c-----------------------------------------------------------------------
28 implicit none
29 #include "SIZE.h"
30 #include "fizhi_SIZE.h"
31 #include "chronos.h"
32
33 C Argument list declarations
34 integer myid,im1,im2,jm1,jm2,idim1,idim2,jdim1,jdim2
35 integer Nrphin,Nsxin,Nsyin,bi,bj,nchp
36 integer nchptot(Nsxin,Nsyin),nchpland(Nsxin,Nsyin)
37 _RL uphy(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
38 _RL vphy(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
39 _RL thphy(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
40 _RL sphy(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
41 _RL pephy(idim1:idim2,jdim1:jdim2,Nrphin+1,Nsxin,Nsyin)
42 _RS lons(idim1:idim2,jdim1:jdim2,Nsxin,Nsyin)
43 _RS lats(idim1:idim2,jdim1:jdim2,Nsxin,Nsyin)
44 _RS phi0surf(idim1:idim2,jdim1:jdim2,Nsxin,Nsyin)
45 _RL ctmt(nchp,Nsxin,Nsyin),xxmt(nchp,Nsxin,Nsyin)
46 _RL yymt(nchp,Nsxin,Nsyin)
47 _RL zetamt(nchp,Nsxin,Nsyin)
48 _RL xlmt(nchp,Nrphin,Nsxin,Nsyin),khmt(nchp,Nrphin,Nsxin,Nsyin)
49 _RL tke(nchp,Nrphin,Nsxin,Nsyin)
50 _RL tgz(im2,jm2,Nsxin,Nsyin)
51 _RL sst(idim1:idim2,jdim1:jdim2,Nsxin,Nsyin)
52 _RL sice(idim1:idim2,jdim1:jdim2,Nsxin,Nsyin)
53 _RL phis_var(im2,jm2,Nsxin,Nsyin)
54 integer landtype(im2,jm2,Nsxin,Nsyin)
55 _RL fracland(im2,jm2,Nsxin,Nsyin),emiss(im2,jm2,10,Nsxin,Nsyin)
56 _RL albvisdr(im2,jm2,Nsxin,Nsyin),albvisdf(im2,jm2,Nsxin,Nsyin)
57 _RL albnirdr(im2,jm2,Nsxin,Nsyin),albnirdf(im2,jm2,Nsxin,Nsyin)
58 _RL chfr(nchp,Nsxin,Nsyin),alai(nchp,Nsxin,Nsyin)
59 _RL agrn(nchp,Nsxin,Nsyin)
60 integer ityp(nchp,Nsxin,Nsyin),igrd(nchp,Nsxin,Nsyin)
61 _RL chlat(nchp,Nsxin,Nsyin),chlon(nchp,Nsxin,Nsyin)
62 _RL tcanopy(nchp,Nsxin,Nsyin),tdeep(nchp,Nsxin,Nsyin)
63 _RL ecanopy(nchp,Nsxin,Nsyin),swetshal(nchp,Nsxin,Nsyin)
64 _RL swetroot(nchp,Nsxin,Nsyin),swetdeep(nchp,Nsxin,Nsyin)
65 _RL snodep(nchp,Nsxin,Nsyin),capac(nchp,Nsxin,Nsyin)
66 _RL o3(im2,jm2,Nrphin,Nsxin,Nsyin)
67 _RL qstr(im2,jm2,Nrphin,Nsxin,Nsyin)
68 _RL co2,cfc11,cfc12,cfc22,n2o(Nrphin),methane(Nrphin)
69
70 integer iras(Nsxin,Nsyin)
71 integer nlwcld(Nsxin,Nsyin),nlwlz(Nsxin,Nsyin)
72 integer nswcld(Nsxin,Nsyin),nswlz(Nsxin,Nsyin)
73 integer imstturbsw(Nsxin,Nsyin),imstturblw(Nsxin,Nsyin)
74 _RL cldtotlwin(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
75 _RL cldraslwin(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
76 _RL cldlsplwin(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
77 _RL lwlzin(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
78 _RL cldtotswin(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
79 _RL cldrasswin(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
80 _RL cldlspswin(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
81 _RL swlzin(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
82 _RL qliqaveswin(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
83 _RL qliqavelwin(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
84 _RL fccaveswin(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
85 _RL fccavelwin(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
86 _RL rainlspin(idim1:idim2,jdim1:jdim2,Nsxin,Nsyin)
87 _RL rainconin(idim1:idim2,jdim1:jdim2,Nsxin,Nsyin)
88 _RL snowfallin(idim1:idim2,jdim1:jdim2,Nsxin,Nsyin)
89
90
91 _RL duphy(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
92 _RL dvphy(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
93 _RL dthphy(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
94 _RL dsphy(idim1:idim2,jdim1:jdim2,Nrphin,Nsxin,Nsyin)
95
96
97 c Local Variables
98 integer ptracer,ntracer
99 parameter (ptracer = 1)
100 parameter (ntracer = 1)
101
102 _RL xlats(sNx,sNy),xlons(sNx,sNy),sea_ice(sNx,sNy)
103 _RL p(sNx,sNy,Nsx,Nsy)
104 _RL u(sNx,sNy,Nrphys),v(sNx,sNy,Nrphys),t(sNx,sNy,Nrphys)
105 _RL q(sNx,sNy,Nrphys,ntracer)
106 _RL pl(sNx,sNy,Nrphys,Nsx,Nsy),pkl(sNx,sNy,Nrphys,Nsx,Nsy)
107 _RL ple(sNx,sNy,Nrphys+1,Nsx,Nsy)
108 _RL pkle(sNx,sNy,Nrphys+1,Nsx,Nsy)
109 _RL dpres(sNx,sNy,Nrphys,Nsx,Nsy)
110 _RL lwdt(sNx,sNy,Nrphys,Nsx,Nsy)
111 _RL lwdtclr(sNx,sNy,Nrphys,Nsx,Nsy)
112 _RL swdt(sNx,sNy,Nrphys,Nsx,Nsy)
113 _RL swdtclr(sNx,sNy,Nrphys,Nsx,Nsy)
114 _RL turbu(sNx,sNy,Nrphys,Nsx,Nsy)
115 _RL turbv(sNx,sNy,Nrphys,Nsx,Nsy)
116 _RL turbt(sNx,sNy,Nrphys,Nsx,Nsy)
117 _RL turbq(sNx,sNy,Nrphys,ntracer,Nsx,Nsy)
118 _RL moistu(sNx,sNy,Nrphys,Nsx,Nsy)
119 _RL moistv(sNx,sNy,Nrphys,Nsx,Nsy)
120 _RL moistt(sNx,sNy,Nrphys,Nsx,Nsy)
121 _RL moistq(sNx,sNy,Nrphys,ntracer,Nsx,Nsy)
122 _RL radswt(sNx,sNy,Nsx,Nsy),radswg(sNx,sNy,Nsx,Nsy)
123 _RL swgclr(sNx,sNy,Nsx,Nsy)
124 _RL fdirpar(sNx,sNy,Nsx,Nsy),fdifpar(sNx,sNy,Nsx,Nsy)
125 _RL osr(sNx,sNy,Nsx,Nsy),osrclr(sNx,sNy,Nsx,Nsy)
126 _RL tg0(sNx,sNy,Nsx,Nsy),radlwg(sNx,sNy,Nsx,Nsy)
127 _RL lwgclr(sNx,sNy,Nsx,Nsy),st4(sNx,sNy,Nsx,Nsy)
128 _RL dst4(sNx,sNy,Nsx,Nsy),dlwdtg(sNx,sNy,Nrphys,Nsx,Nsy)
129 _RL qq(sNx,sNy,Nrphys,Nsx,Nsy)
130 integer i,j,L
131 _RL getcon, kappa, p0kappa, s0, ra
132 _RL cosz(sNx,sNy)
133 _RL cldtot_lw(sNx,sNy,Nrphys)
134 _RL cldras_lw(sNx,sNy,Nrphys)
135 _RL cldlsp_lw(sNx,sNy,Nrphys)
136 _RL lwlz(sNx,sNy,Nrphys)
137 _RL cldtot_sw(sNx,sNy,Nrphys)
138 _RL cldras_sw(sNx,sNy,Nrphys)
139 _RL cldlsp_sw(sNx,sNy,Nrphys)
140 _RL swlz(sNx,sNy,Nrphys)
141 _RL qliqavesw(sNx,sNy,Nrphys)
142 _RL qliqavelw(sNx,sNy,Nrphys)
143 _RL fccavesw(sNx,sNy,Nrphys)
144 _RL fccavelw(sNx,sNy,Nrphys)
145 _RL rainlsp(sNx,sNy)
146 _RL raincon(sNx,sNy)
147 _RL snowfall(sNx,sNy)
148
149 _RL tempij(sNx,sNy)
150 _RL tempi(2)
151
152 _RL kF,sigma_b,ks,ka,deg2rad,pi,atm_po,atm_kappa,termp,kv,kT
153 _RL term1,term2,thetalim,thetaeq,recip_p0g
154
155 logical alarm
156 external alarm
157
158 C***********************************************************************
159
160 kF=1. _d 0/86400. _d 0
161 sigma_b = 0.7 _d 0
162 ka=1. _d 0/(40. _d 0*86400. _d 0)
163 ks=1. _d 0/(4. _d 0 *86400. _d 0)
164 pi = getcon('PI')
165 atm_kappa = getcon('KAPPA')
166 atm_po = getcon('ATMPOPA')
167 deg2rad = getcon('DEG2RAD')
168
169 do L = 1,Nrphys
170 do j = jm1,jm2
171 do i = im1,im2
172 recip_P0g= 1. _d 0 / pephy(i,j,Nrphys+1,bi,bj)
173 c U and V terms:
174 termP=0.5 _d 0*((pephy(i,j,L,bi,bj)+pephy(i,j,L+1,bi,bj))
175 & *recip_P0g )
176 kV=kF*MAX( 0. _d 0, (termP-sigma_b)/(1. _d 0-sigma_b) )
177 duphy(i,j,L,bi,bj)= -kV*uphy(i,j,L,bi,bj)
178 dvphy(i,j,L,bi,bj)= -kV*vphy(i,j,L,bi,bj)
179
180 c T terms
181 C-- Forcing term(s)
182 term1=60. _d 0*(sin(lats(I,J,bi,bj)*deg2rad)**2)
183 termP=0.5 _d 0*( pephy(i,j,L,bi,bj) + pephy(i,j,L+1,bi,bj) )
184 term2=10. _d 0*log(termP/atm_po)
185 & *(cos(lats(I,J,bi,bj)*deg2rad)**2)
186 thetaLim = 200. _d 0/ ((termP/atm_po)**atm_kappa)
187 thetaEq=315. _d 0-term1-term2
188 thetaEq=MAX(thetaLim,thetaEq)
189 kT=ka+(ks-ka)
190 & *MAX(0. _d 0,
191 & (termP*recip_P0g-sigma_b)/(1. _d 0-sigma_b) )
192 & *COS((lats(I,J,bi,bj)*deg2rad))**4
193 if(termP*recip_P0g.gt.0.04)then
194 dthphy(i,j,L,bi,bj)=- kT*( thphy(I,J,L,bi,bj)-thetaEq )
195 else
196 dthphy(i,j,L,bi,bj)=0.
197 endif
198
199 c S terms (hs runs dry - no moisture)
200 C-- Forcing term(s)
201 dsphy(i,j,L,bi,bj)=0.
202
203 enddo
204 enddo
205 enddo
206
207 return
208 end

  ViewVC Help
Powered by ViewVC 1.1.22