/[MITgcm]/MITgcm/model/src/ini_p_ground.F
ViewVC logotype

Contents of /MITgcm/model/src/ini_p_ground.F

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


Revision 1.1 - (show annotations) (download)
Fri Jul 6 21:48:31 2001 UTC (22 years, 10 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint40pre3, checkpoint40pre7, checkpoint40pre6, checkpoint40pre9, checkpoint40pre8, checkpoint40pre2, checkpoint40pre4, checkpoint40pre5, checkpoint40
read orography (in m) from file and convert to Surface Pressure (Ro_surf)

1 C $Header: $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5
6 SUBROUTINE INI_P_GROUND(
7 I Hfld,
8 O Pfld,
9 I myThid )
10 C /==========================================================\
11 C | SUBROUTINE INI_P_GROUND |
12 C | o Convert Topography [m] to (reference) Surface Pressure |
13 C | according to tRef profile, |
14 C | using same discretisation as in calc_phi_hyd |
15 C | |
16 C \==========================================================/
17 IMPLICIT NONE
18 C == Global variables ==
19 #include "SIZE.h"
20 #include "GRID.h"
21 #include "EEPARAMS.h"
22 #include "PARAMS.h"
23
24 C == Routine arguments ==
25 C Hfld (input) : Ground elevation [m]
26 C Pfld (output) : reference Pressure at the ground [Pa]
27 _RS Hfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
28 _RS Pfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
29 INTEGER myThid
30
31 c #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
32
33 C == Local variables ==
34 INTEGER bi,bj,i,j,K, Ks
35 _RL ddPI, Po_surf
36 _RL phiRef(2*Nr+1), rHalf(2*Nr+1)
37
38 DO K=1,Nr
39 rHalf(2*K-1) = rF(K)
40 rHalf(2*K) = rC(K)
41 ENDDO
42 rHalf(2*Nr+1) = rF(Nr+1)
43
44 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
45 C- Compute Reference Geopotential at Half levels :
46 C Tracer level: phiRef(2K) ; Interface_W level: phiRef(2K+1)
47
48 phiRef(1) = 0.
49
50 IF (Integr_GeoPot.EQ.1) THEN
51 C- Finite Volume Form, linear by half level :
52 DO K=1,2*Nr
53 Ks = (K+1)/2
54 ddPI=atm_cp*( ((rHalf( K )/atm_po)**atm_kappa)
55 & -((rHalf(K+1)/atm_po)**atm_kappa) )
56 phiRef(K+1) = phiRef(K)+ddPI*tRef(Ks)
57 ENDDO
58 C------
59 ELSE
60 C- Finite Difference Form, linear between Tracer level :
61 C works with Integr_GeoPot = 0, 2 or 3
62 K = 1
63 ddPI=atm_cp*( ((rF(K)/atm_po)**atm_kappa)
64 & -((rC(K)/atm_po)**atm_kappa) )
65 phiRef(2*K) = phiRef(1) + ddPI*tRef(K)
66 DO K=1,Nr-1
67 ddPI=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
68 & -((rC(K+1)/atm_po)**atm_kappa) )
69 phiRef(2*K+1) = phiRef(2*K) + ddPI*0.5*tRef(K)
70 phiRef(2*K+2) = phiRef(2*K)
71 & + ddPI*0.5*(tRef(K)+tRef(K+1))
72 ENDDO
73 K = Nr
74 ddPI=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
75 & -((rF(K+1)/atm_po)**atm_kappa) )
76 phiRef(2*K+1) = phiRef(2*K) + ddPI*tRef(K)
77 C------
78 ENDIF
79
80 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
81
82 C- Convert phiRef to Z unit :
83 DO K=1,2*Nr+1
84 phiRef(K) = phiRef(K)*recip_gravity
85 ENDDO
86
87 C- Write to check :
88 write(6,'(2A)') 'INI_P_GROUND: PhiRef/g [m] at Center ',
89 & '(integer) and Interface (half-int.) levels:'
90 DO K=1,2*Nr+1
91 write(6,'(A,F5.1,A,F10.1,A,F12.3)') ' K=', K*0.5,
92 & ' ; r=', rHalf(K), ' ; phiRef/g=', phiRef(K)
93 ENDDO
94
95 C- Find Po_surf : Linear between 2 half levels :
96 DO bj = myByLo(myThid), myByHi(myThid)
97 DO bi = myBxLo(myThid), myBxHi(myThid)
98 DO j=1,sNy
99 DO i=1,sNx
100 Ks = 1
101 DO K=1,2*Nr
102 IF (Hfld(i,j,bi,bj).GE.phiRef(K)) Ks = K
103 ENDDO
104 Po_surf = rHalf(Ks) + (rHalf(Ks+1)-rHalf(Ks))*
105 & (Hfld(i,j,bi,bj)-phiRef(Ks))/(phiRef(Ks+1)-phiRef(Ks))
106
107 c IF (Hfld(i,j,bi,bj).LT.phiRef(1)) Po_surf= rHalf(1)
108 c IF (Hfld(i,j,bi,bj).GT.phiRef(2*Nr+1)) Po_surf= rHalf(2*Nr+1)
109 Pfld(i,j,bi,bj) = Po_surf
110 ENDDO
111 ENDDO
112 ENDDO
113 ENDDO
114
115 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
116
117 c #endif/* INCLUDE_PHIHYD_CALCULATION_CODE */
118
119 RETURN
120 END

  ViewVC Help
Powered by ViewVC 1.1.22