pkg/aim/phy_shtorh.F

C $Header: /u/gcmpack/models/MITgcmUV/pkg/aim/phy_shtorh.F,v 1.3 2001/05/29 19:28:53 cnh Exp $
C $Name:  $
 
      SUBROUTINE SHTORH (IMODE,NGP,TA,PS,SIG,QA,RH,QSAT,myThid)
C--
C--   SUBROUTINE SHTORH (IMODE,NGP,TA,PS,SIG,QA,RH,QSAT)
C--
C--   Purpose: compute saturation specific humidity and 
C--            relative hum. from specific hum. (or viceversa)
C--   Input:   IMODE  : mode of operation
C--            NGP    : no. of grid-points
C--            TA     : abs. temperature
C--            PS     : normalized pressure   (=  p/1000_hPa) [if SIG < 0]
C--                   : normalized sfc. pres. (= ps/1000_hPa) [if SIG > 0]
C--            SIG    : sigma level
C--            QA     : specific humidity in g/kg [if IMODE > 0]
C--            RH     : relative humidity         [if IMODE < 0]
C--            QSAT   : saturation spec. hum. in g/kg 
C--   Output:  RH     : relative humidity         [if IMODE > 0] 
C--            QA     : specific humidity in g/kg [if IMODE < 0]
C--        


      IMPLICIT rEAL*8 (A-H,O-Z)
      INTEGER IMODE,NGP,J
      INTEGER  myThid


CcnhDebugStarts
#include "SIZE.h"
CcnhDebugEnds
      REAL TA(NGP), PS(NGP), QA(NGP), RH(NGP), QSAT(NGP)
C
C---  1. Compute Qsat (g/kg) from T (degK) and normalized pres. P (= p/1000_hPa)
C        If SIG > 0, P = Ps * sigma, otherwise P = Ps(1) = const. 
C
      E0=  6.108 _d -3
      C1= 17.269 _d 0
      C2= 21.875 _d 0
      T0=273.16 _d 0
      T1= 35.86 _d 0
      T2=  7.66 _d 0
C 
      DO 110 J=1,NGP
        QSAT(J)=0.
        IF (TA(J).GE.T0) THEN
          QSAT(J)=E0*EXP(C1*(TA(J)-T0)/(TA(J)-T1))
        ELSE IF ( TA(J).GT.0.) then
          QSAT(J)=E0*EXP(C2*(TA(J)-T0)/(TA(J)-T2))
        ENDIF
  110 CONTINUE
C
      IF (SIG.LE.0.0) THEN
        DO 120 J=1,NGP
          QSAT(J)=622. _d 0*QSAT(J)/(PS(1)-0.378 _d 0*QSAT(J))
  120   CONTINUE
      ELSE
        DO 130 J=1,NGP
          QSAT(J)=622. _d 0*QSAT(J)/(SIG*PS(J)-0.378 _d 0*QSAT(J))
  130   CONTINUE
      ENDIF
chh      write(0,*) 'MAXVAL(QSAT)=',MAXVAL(QSAT)
chh      write(0,*) 'MINVAL(QSAT)=',MINVAL(QSAT)
C
C---  2. Compute rel.hum. RH=Q/Qsat (IMODE>0), or Q=RH*Qsat (IMODE<0)
C
      IF (IMODE.GT.0) THEN
        DO 210 J=1,NGP
          IF(QSAT(J).ne.0.) then
            RH(J)=QA(J)/QSAT(J)
          ELSE
            RH(J)=0.
          ENDIF
  210   CONTINUE
      ELSE IF (IMODE.LT.0) THEN
        DO 220 J=1,NGP
          QA(J)=RH(J)*QSAT(J)
  220   CONTINUE
      ENDIF
chh      write(0,*) 'MAXVAL(QA)=',MAXVAL(QA)
chh      write(0,*) 'MINVAL(QA)=',MINVAL(QA)
chh      write(0,*) 'MAXVAL(RH)=',MAXVAL(RH)
chh      write(0,*) 'MINVAL(RH)=',MINVAL(RH)
C                               
      RETURN
      END

      SUBROUTINE ZMEDDY (NLON,NLAT,FF,ZM,EDDY)


      IMPLICIT rEAL*8 (A-H,O-Z)
      INTEGER NLON,NLAT,I,J


C
C *** Decompose a field into zonal-mean and eddy component
C
      REAL FF(NLON,NLAT), ZM(NLAT), EDDY(NLON,NLAT)
C
      RNLON=1./NLON
C
      DO 130 J=1,NLAT
C
        ZM(J)=0.
        DO 110 I=1,NLON
          ZM(J)=ZM(J)+FF(I,J)
 110    CONTINUE
        ZM(J)=ZM(J)*RNLON
C
        DO 120 I=1,NLON
          EDDY(I,J)=FF(I,J)-ZM(J)
 120    CONTINUE
C
 130  CONTINUE
C
C--
      RETURN
      END
C
1	C $Header: /u/gcmpack/models/MITgcmUV/pkg/aim/phy_shtorh.F,v 1.3 2001/05/29 19:28:53 cnh Exp $
2	C $Name: $
3
4	SUBROUTINE SHTORH (IMODE,NGP,TA,PS,SIG,QA,RH,QSAT,myThid)
5	C--
6	C-- SUBROUTINE SHTORH (IMODE,NGP,TA,PS,SIG,QA,RH,QSAT)
7	C--
8	C-- Purpose: compute saturation specific humidity and
9	C-- relative hum. from specific hum. (or viceversa)
10	C-- Input: IMODE : mode of operation
11	C-- NGP : no. of grid-points
12	C-- TA : abs. temperature
13	C-- PS : normalized pressure (= p/1000_hPa) [if SIG < 0]
14	C-- : normalized sfc. pres. (= ps/1000_hPa) [if SIG > 0]
15	C-- SIG : sigma level
16	C-- QA : specific humidity in g/kg [if IMODE > 0]
17	C-- RH : relative humidity [if IMODE < 0]
18	C-- QSAT : saturation spec. hum. in g/kg
19	C-- Output: RH : relative humidity [if IMODE > 0]
20	C-- QA : specific humidity in g/kg [if IMODE < 0]
21	C--
22
23
24	IMPLICIT rEAL*8 (A-H,O-Z)
25	INTEGER IMODE,NGP,J
26	INTEGER myThid
27
28
29	CcnhDebugStarts
30	#include "SIZE.h"
31	CcnhDebugEnds
32	REAL TA(NGP), PS(NGP), QA(NGP), RH(NGP), QSAT(NGP)
33	C
34	C--- 1. Compute Qsat (g/kg) from T (degK) and normalized pres. P (= p/1000_hPa)
35	C If SIG > 0, P = Ps * sigma, otherwise P = Ps(1) = const.
36	C
37	E0= 6.108 _d -3
38	C1= 17.269 _d 0
39	C2= 21.875 _d 0
40	T0=273.16 _d 0
41	T1= 35.86 _d 0
42	T2= 7.66 _d 0
43	C
44	DO 110 J=1,NGP
45	QSAT(J)=0.
46	IF (TA(J).GE.T0) THEN
47	QSAT(J)=E0EXP(C1(TA(J)-T0)/(TA(J)-T1))
48	ELSE IF ( TA(J).GT.0.) then
49	QSAT(J)=E0EXP(C2(TA(J)-T0)/(TA(J)-T2))
50	ENDIF
51	110 CONTINUE
52	C
53	IF (SIG.LE.0.0) THEN
54	DO 120 J=1,NGP
55	QSAT(J)=622. _d 0QSAT(J)/(PS(1)-0.378 _d 0QSAT(J))
56	120 CONTINUE
57	ELSE
58	DO 130 J=1,NGP
59	QSAT(J)=622. _d 0QSAT(J)/(SIGPS(J)-0.378 _d 0*QSAT(J))
60	130 CONTINUE
61	ENDIF
62	chh write(0,*) 'MAXVAL(QSAT)=',MAXVAL(QSAT)
63	chh write(0,*) 'MINVAL(QSAT)=',MINVAL(QSAT)
64	C
65	C--- 2. Compute rel.hum. RH=Q/Qsat (IMODE>0), or Q=RH*Qsat (IMODE<0)
66	C
67	IF (IMODE.GT.0) THEN
68	DO 210 J=1,NGP
69	IF(QSAT(J).ne.0.) then
70	RH(J)=QA(J)/QSAT(J)
71	ELSE
72	RH(J)=0.
73	ENDIF
74	210 CONTINUE
75	ELSE IF (IMODE.LT.0) THEN
76	DO 220 J=1,NGP
77	QA(J)=RH(J)*QSAT(J)
78	220 CONTINUE
79	ENDIF
80	chh write(0,*) 'MAXVAL(QA)=',MAXVAL(QA)
81	chh write(0,*) 'MINVAL(QA)=',MINVAL(QA)
82	chh write(0,*) 'MAXVAL(RH)=',MAXVAL(RH)
83	chh write(0,*) 'MINVAL(RH)=',MINVAL(RH)
84	C
85	RETURN
86	END
87
88	SUBROUTINE ZMEDDY (NLON,NLAT,FF,ZM,EDDY)
89
90
91	IMPLICIT rEAL*8 (A-H,O-Z)
92	INTEGER NLON,NLAT,I,J
93
94
95	C
96	C *** Decompose a field into zonal-mean and eddy component
97	C
98	REAL FF(NLON,NLAT), ZM(NLAT), EDDY(NLON,NLAT)
99	C
100	RNLON=1./NLON
101	C
102	DO 130 J=1,NLAT
103	C
104	ZM(J)=0.
105	DO 110 I=1,NLON
106	ZM(J)=ZM(J)+FF(I,J)
107	110 CONTINUE
108	ZM(J)=ZM(J)*RNLON
109	C
110	DO 120 I=1,NLON
111	EDDY(I,J)=FF(I,J)-ZM(J)
112	120 CONTINUE
113	C
114	130 CONTINUE
115	C
116	C--
117	RETURN
118	END
119	C