1 |
subroutine do_fizhi(uphy,vphy,thphy,sphy,pephy, |
2 |
. ctmt,xxmt,yymt,zetamt,xlmt,khmt,tke, |
3 |
. xC,yC, |
4 |
.im1,im2,jm1,jm2,Nrphys,Nsx,Nsy,idim1,idim2,jdim1,jdim2,bi,bj,nchp, |
5 |
. duphy,dvphy,dthphy,dsphy) |
6 |
c----------------------------------------------------------------------- |
7 |
c Dummy routine to calculate physics increments - here set them to |
8 |
c the Held-Suarez forcing terms |
9 |
c----------------------------------------------------------------------- |
10 |
implicit none |
11 |
#include "CPP_OPTIONS.h" |
12 |
|
13 |
integer im1,im2,jm1,jm2,idim1,idim2,jdim1,jdim2 |
14 |
integer Nrphys,Nsx,Nsy,bi,bj,nchp |
15 |
_RL uphy(im1:im2,jm1:jm2,Nrphys,Nsx,Nsy) |
16 |
_RL vphy(im1:im2,jm1:jm2,Nrphys,Nsx,Nsy) |
17 |
_RL thphy(im1:im2,jm1:jm2,Nrphys,Nsx,Nsy) |
18 |
_RL sphy(im1:im2,jm1:jm2,Nrphys,Nsx,Nsy) |
19 |
_RL pephy(im1:im2,jm1:jm2,Nrphys+1,Nsx,Nsy) |
20 |
_RL ctmt(nchp,Nsx,Nsy),xxmt(nchp,Nsx,Nsy),yymt(nchp,Nsx,Nsy) |
21 |
_RL zetamt(nchp,Nsx,Nsy) |
22 |
_RL xlmt(nchp,Nrphys,Nsx,Nsy),khmt(nchp,Nrphys,Nsx,Nsy) |
23 |
_RL tke(nchp,Nrphys,Nsx,Nsy) |
24 |
_RL duphy(im1:im2,jm1:jm2,Nrphys,Nsx,Nsy) |
25 |
_RL dvphy(im1:im2,jm1:jm2,Nrphys,Nsx,Nsy) |
26 |
_RL dthphy(im1:im2,jm1:jm2,Nrphys,Nsx,Nsy) |
27 |
_RL dsphy(im1:im2,jm1:jm2,Nrphys,Nsx,Nsy) |
28 |
_RS xC(im1:im2,jm1:jm2,Nsx,Nsy) |
29 |
_RS yC(im1:im2,jm1:jm2,Nsx,Nsy) |
30 |
c |
31 |
integer i,j,L |
32 |
_RL kF,sigma_b,ks,ka,deg2rad,pi,atm_po,atm_kappa,termp,kv,kT |
33 |
_RL term1,term2,thetalim,thetaeq,recip_p0g |
34 |
_RL getcon |
35 |
|
36 |
kF=1. _d 0/86400. _d 0 |
37 |
sigma_b = 0.7 _d 0 |
38 |
ka=1. _d 0/(40. _d 0*86400. _d 0) |
39 |
ks=1. _d 0/(4. _d 0 *86400. _d 0) |
40 |
pi = getcon('PI') |
41 |
atm_kappa = getcon('KAPPA') |
42 |
atm_po = getcon('ATMPOPA') |
43 |
deg2rad = getcon('DEG2RAD') |
44 |
|
45 |
do L = 1,Nrphys |
46 |
do j = jdim1,jdim2 |
47 |
do i = idim1,idim2 |
48 |
recip_P0g= 1. _d 0 / pephy(i,j,1,bi,bj) |
49 |
c U and V terms: |
50 |
termP=0.5 _d 0*((pephy(i,j,L,bi,bj)+pephy(i,j,L+1,bi,bj)) |
51 |
& *recip_P0g ) |
52 |
kV=kF*MAX( 0. _d 0, (termP-sigma_b)/(1. _d 0-sigma_b) ) |
53 |
duphy(i,j,L,bi,bj)= -kV*uphy(i,j,L,bi,bj) |
54 |
dvphy(i,j,L,bi,bj)= -kV*vphy(i,j,L,bi,bj) |
55 |
|
56 |
c T terms |
57 |
C-- Forcing term(s) |
58 |
term1=60. _d 0*(sin(yC(I,J,bi,bj)*deg2rad)**2) |
59 |
termP=0.5 _d 0*( pephy(i,j,L,bi,bj) + pephy(i,j,L+1,bi,bj) ) |
60 |
term2=10. _d 0*log(termP/atm_po) |
61 |
& *(cos(yC(I,J,bi,bj)*deg2rad)**2) |
62 |
thetaLim = 200. _d 0/ ((termP/atm_po)**atm_kappa) |
63 |
thetaEq=315. _d 0-term1-term2 |
64 |
thetaEq=MAX(thetaLim,thetaEq) |
65 |
kT=ka+(ks-ka) |
66 |
& *MAX(0. _d 0, |
67 |
& (termP*recip_P0g-sigma_b)/(1. _d 0-sigma_b) ) |
68 |
& *COS((yC(I,J,bi,bj)*deg2rad))**4 |
69 |
if(termP*recip_P0g.gt.0.04)then |
70 |
dthphy(i,j,L,bi,bj)=- kT*( thphy(I,J,L,bi,bj)-thetaEq ) |
71 |
else |
72 |
dthphy(i,j,L,bi,bj)=0. |
73 |
endif |
74 |
|
75 |
c S terms (hs runs dry - no moisture) |
76 |
C-- Forcing term(s) |
77 |
dsphy(i,j,L,bi,bj)=0. |
78 |
|
79 |
enddo |
80 |
enddo |
81 |
enddo |
82 |
|
83 |
return |
84 |
end |