/[MITgcm]/MITgcm/pkg/aim/phy_suflux.F
ViewVC logotype

Diff of /MITgcm/pkg/aim/phy_suflux.F

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

revision 1.2 by adcroft, Fri Feb 2 21:36:29 2001 UTC revision 1.3 by cnh, Tue May 29 19:28:53 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4        SUBROUTINE SUFLUX (PSA,UA,VA,TA,QA,RH,QSAT,PHI,        SUBROUTINE SUFLUX (PSA,UA,VA,TA,QA,RH,QSAT,Vsurfsq,PHI,
5       *                   PHI0,FMASK,TLAND,TSEA,SWAV,SSR,SLR,       *                   PHI0,FMASK,TLAND,TSEA,SWAV,SSR,SLR,
6       *                   USTR,VSTR,SHF,EVAP,T0,Q0,QSAT0,SPEED0)       *                   DRAG,USTR,VSTR,SHF,EVAP,T0,Q0,QSAT0,SPEED0,
7         &                   myThid)
8  C--  C--
9  C--   SUBROUTINE SUFLUX (PSA,UA,VA,TA,QA,RH,PHI,  C--   SUBROUTINE SUFLUX (PSA,UA,VA,TA,QA,RH,PHI,
10  C--  *                   PHI0,FMASK,TLAND,TSEA,SWAV,  C--  *                   PHI0,FMASK,TLAND,TSEA,SWAV,
# Line 26  C--   Output:  USTR   = u stress Line 27  C--   Output:  USTR   = u stress
27  C--            VSTR   = v stress                         (2-dim)  C--            VSTR   = v stress                         (2-dim)
28  C--            SHF    = sensible heat flux               (2-dim)  C--            SHF    = sensible heat flux               (2-dim)
29  C--            EVAP   = evaporation [g/(m^2 s)]          (2-dim)  C--            EVAP   = evaporation [g/(m^2 s)]          (2-dim)
30    C - added (jmc) :
31    C--           Vsurfsq = square of surface wind speed (2-dim,input)
32    C--                       ==> UA,VA are no longer used
33    C--            DRAG   = surface Drag term (= Cd*|V|) (2-dim,output)
34    C--                       ==> USTR,VSTR are no longer used
35    C--         myThid    = Instance number of this instance of the
36    C--                     the routine.
37  C--  C--
38    
39    
# Line 33  C-- Line 41  C--
41    
42    
43  C     Resolution parameters  C     Resolution parameters
44    #include "EEPARAMS.h"
45  #include "atparam.h"  #include "atparam.h"
46  #include "atparam1.h"  #include "atparam1.h"
47  #include "Lev_def.h"  #include "Lev_def.h"
# Line 46  C     Physical constants + functions of Line 55  C     Physical constants + functions of
55  C     Surface flux constants  C     Surface flux constants
56    
57  #include "com_sflcon.h"  #include "com_sflcon.h"
58    
59  C  C
60        REAL PSA(NGP), UA(NGP,NLEV), VA(NGP,NLEV), TA(NGP,NLEV),        REAL PSA(NGP), UA(NGP,NLEV), VA(NGP,NLEV), TA(NGP,NLEV),
61       *     QA(NGP,NLEV), RH(NGP,NLEV), QSAT(NGP,NLEV), PHI(NGP,NLEV),       *     QA(NGP,NLEV), RH(NGP,NLEV), QSAT(NGP,NLEV), PHI(NGP,NLEV),
62       *     PHI0(NGP), FMASK(NGP), TLAND(NGP), TSEA(NGP), SWAV(NGP),       *     PHI0(NGP), FMASK(NGP), TLAND(NGP), TSEA(NGP), SWAV(NGP),
63       *     SSR(NGP), SLR(NGP)       *     SSR(NGP), SLR(NGP)
64          REAL Vsurfsq(NGP), DRAG(NGP)
65          INTEGER myThid
66  C  C
67        REAL USTR(NGP,3), VSTR(NGP,3), SHF(NGP,3), EVAP(NGP,3)        REAL USTR(NGP,3), VSTR(NGP,3), SHF(NGP,3), EVAP(NGP,3)
68  C                                                                        C                                                                      
# Line 79  Cchdbg Line 91  Cchdbg
91        SAVE denmoy1, denmoy2        SAVE denmoy1, denmoy2
92        REAL Q0moy(NGP)        REAL Q0moy(NGP)
93        SAVE Q0moy        SAVE Q0moy
94          REAL factwind2
95  C  C
96  C--   1. Extrapolation of wind, temp, hum. and density to the surface  C--   1. Extrapolation of wind, temp, hum. and density to the surface
97  C  C
98  C     1.1 Wind components  C     1.1 Wind components
99  C    C  
100        DO J=1,NGP        DO J=1,NGP
101          IF ( NLEVxyU(J) .GT. 0 ) THEN  c_jmc   IF ( NLEVxyU(J) .GT. 0 ) THEN
102           U0(J) = FWIND0*UA(J,NLEVxyU(J))  c_jmc    U0(J) = FWIND0*UA(J,NLEVxyU(J))
103          ENDIF  c_jmc   ENDIF
104          IF ( NLEVxyV(J) .GT. 0 ) THEN  c_jmc   IF ( NLEVxyV(J) .GT. 0 ) THEN
105           V0(J) = FWIND0*VA(J,NLEVxyV(J))  c_jmc    V0(J) = FWIND0*VA(J,NLEVxyV(J))
106    c_jmc   ENDIF
107    c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
108    C-- Note(jmc): If we interpolate U,V from C_grid to A_grid (center of the mesh)
109    C     the above distinction between NLEVxy_Upt,Vpt is no longer valid !
110    C     same thing (again) few lines below.
111    c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
112            IF ( NLEVxy(J,myThid) .GT. 0 ) THEN
113             U0(J) = FWIND0*UA(J,NLEVxy(J,myThid))
114             V0(J) = FWIND0*VA(J,NLEVxy(J,myThid))
115            ELSE
116             U0(J) = 0.0
117             V0(J) = 0.0
118          ENDIF          ENDIF
119  C       U0(J) = UA(J,5)  C       U0(J) = UA(J,5)
120  C       V0(J) = VA(J,5)  C       V0(J) = VA(J,5)
# Line 100  C Line 125  C
125        GTEMP0 = 1.D0-FTEMP0        GTEMP0 = 1.D0-FTEMP0
126        RCP = 1.D0/CP        RCP = 1.D0/CP
127        DO J=1,NGP        DO J=1,NGP
128          NL1(J)=NLEVxy(J)-1          NL1(J)=NLEVxy(J,myThid)-1
129        ENDDO        ENDDO
130  C  C
131        DO J=1,NGP        DO J=1,NGP
132         IF ( NLEVxy(J) .GT. 0 ) THEN         IF ( NLEVxy(J,myThid) .GT. 0 ) THEN
133          T0(J,1) = TA(J,NLEVxy(J))+WVI(NLEVxy(J),2)*          T0(J,1) = TA(J,NLEVxy(J,myThid))+WVI(NLEVxy(J,myThid),2)*
134       &                          (TA(J,NLEVxy(J))-TA(J,NL1(J)))       &                          (TA(J,NLEVxy(J,myThid))-TA(J,NL1(J)))
135  ccccc        T0(J,2) = TA(J,NLEVxy(J))+RCP*(PHI(J,NLEVxy(J))-PHI0(J))  ccccc        T0(J,2) = TA(J,NLEVxy(J))+RCP*(PHI(J,NLEVxy(J))-PHI0(J))
136          T0(J,2) = TA(J,NLEVxy(J))*          T0(J,2) = TA(J,NLEVxy(J,myThid))*
137       &         ((Psurfw(NLEVxy(J))/Psurfs(Nlevxy(J)))**(RD/CP))       &         ((Psurfw(NLEVxy(J,myThid))/
138         &           Psurfs(Nlevxy(J,myThid)))**(RD/CP))
139         ELSE         ELSE
140          T0(J,1) = 273.16 _d 0          T0(J,1) = 273.16 _d 0
141          T0(J,2) = 273.16 _d 0          T0(J,2) = 273.16 _d 0
# Line 125  C Line 151  C
151        GHUM0 = 1. _d 0-FHUM0        GHUM0 = 1. _d 0-FHUM0
152  C  C
153        DO J=1,NGP        DO J=1,NGP
154         IF ( NLEVxy(J) .GT. 0 ) THEN         IF ( NLEVxy(J,myThid) .GT. 0 ) THEN
155          WORK(J)=RH(J,Nlevxy(J))          WORK(J)=RH(J,Nlevxy(J,myThid))
156  cchdbg  cchdbg
157  c        WORK(J) = RH(J,NLEVxy(J))+WVI(NLEVxy(J),2)*  c        WORK(J) = RH(J,NLEVxy(J))+WVI(NLEVxy(J),2)*
158  c     &                          (RH(J,NLEVxy(J))-RH(J,NL1(J)))  c     &                          (RH(J,NLEVxy(J))-RH(J,NL1(J)))
# Line 134  cchdbg Line 160  cchdbg
160         ENDIF         ENDIF
161        ENDDO        ENDDO
162  C  C
163        CALL SHTORH (-1,NGP,T0,PSA,1. _d 0,Q0,WORK,QSAT0)        CALL SHTORH (-1,NGP,T0,PSA,1. _d 0,Q0,WORK,QSAT0,myThid)
164  C  C
165        DO J=1,NGP        DO J=1,NGP
166         IF ( NLEVxy(J) .GT. 0 ) THEN         IF ( NLEVxy(J,myThid) .GT. 0 ) THEN
167          Q0(J)=FHUM0*Q0(J)+GHUM0*QA(J,NLEVxy(J))          Q0(J)=FHUM0*Q0(J)+GHUM0*QA(J,NLEVxy(J,myThid))
168         ENDIF         ENDIF
169        ENDDO        ENDDO
170    
# Line 146  C     1.4 Density * wind speed (includin Line 172  C     1.4 Density * wind speed (includin
172    
173        PRD = P0/RD        PRD = P0/RD
174        VG2 = VGUST*VGUST        VG2 = VGUST*VGUST
175          factwind2 = FWIND0*FWIND0
176  C  C
177        DO J=1,NGP        DO J=1,NGP
178          DENVV(J,3)=(PRD*PSA(J)/T0(J,1))*  c_jmc   SPEED0(J)=SQRT(U0(J)*U0(J)+V0(J)*V0(J)+VG2)
179       *           SQRT(U0(J)*U0(J)+V0(J)*V0(J)+VG2)          SPEED0(J)=SQRT(factwind2*Vsurfsq(J)+VG2)
180            DENVV(J,3)=(PRD*PSA(J)/T0(J,1))*SPEED0(J)
181    c_jmc   DENVV(J,3)=(PRD*PSA(J)/T0(J,1))*
182    c_jmc*           SQRT(U0(J)*U0(J)+V0(J)*V0(J)+VG2)
183  cchdbg        DENVV(J,3)=0.7*(PRD*PSA(J)/T0(J,1))*  cchdbg        DENVV(J,3)=0.7*(PRD*PSA(J)/T0(J,1))*
184  cchdbg     *           SQRT(U0(J)*U0(J)+V0(J)*V0(J)+VG2)  cchdbg     *           SQRT(U0(J)*U0(J)+V0(J)*V0(J)+VG2)
         SPEED0(J)=SQRT(U0(J)*U0(J)+V0(J)*V0(J)+VG2)  
185        ENDDO        ENDDO
186  C  C
187  C     1.5 Stability correction for heat/moisture fluxes  C     1.5 Stability correction for heat/moisture fluxes
# Line 176  C Line 205  C
205  C     2.1 Wind stress  C     2.1 Wind stress
206  C  C
207        DO J=1,NGP        DO J=1,NGP
208         IF ( NLEVxyu(J) .GT. 0  ) THEN  c_jmc  IF ( NLEVxyu(J) .GT. 0  ) THEN
209          USTR(J,1) = -CDL*DENVV(J,3)*UA(J,NLEVxyu(J))  c_jmc   USTR(J,1) = -CDL*DENVV(J,3)*UA(J,NLEVxyu(J))
210          USTR(J,2) = -CDS*DENVV(J,3)*UA(J,NLEVxyu(J))  c_jmc   USTR(J,2) = -CDS*DENVV(J,3)*UA(J,NLEVxyu(J))
211    c_jmc  ENDIF
212    c_jmc  IF ( NLEVxyv(J) .GT. 0  ) THEN
213    c_jmc   VSTR(J,1) = -CDL*DENVV(J,3)*VA(J,NLEVxyv(J))
214    c_jmc   VSTR(J,2) = -CDS*DENVV(J,3)*VA(J,NLEVxyv(j))
215    c_jmc  ENDIF
216    C-- cf Note(jmc) above.
217           IF ( NLEVxy(J,myThid) .GT. 0  ) THEN
218            USTR(J,1) = -CDL*DENVV(J,3)*UA(J,NLEVxy(J,myThid))
219            USTR(J,2) = -CDS*DENVV(J,3)*UA(J,NLEVxy(J,myThid))
220            VSTR(J,1) = -CDL*DENVV(J,3)*VA(J,NLEVxy(J,myThid))
221            VSTR(J,2) = -CDS*DENVV(J,3)*VA(J,NLEVxy(J,myThid))
222           ELSE
223            USTR(J,1) = 0.
224            USTR(J,2) = 0.
225            VSTR(J,1) = 0.
226            VSTR(J,2) = 0.
227         ENDIF         ENDIF
228         IF ( NLEVxyv(J) .GT. 0  ) THEN        ENDDO
229          VSTR(J,1) = -CDL*DENVV(J,3)*VA(J,NLEVxyv(J))  c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
230          VSTR(J,2) = -CDS*DENVV(J,3)*VA(J,NLEVxyv(j))  c   Compute surface drag term (= C_drag*|V| ) to allow direct computation
231    c     of surface stress on C-grid.
232    c     add Land + Sea contributions ; Convert to surface pressure level
233          DO J=1,NGP
234           IF ( NLEVxy(J,myThid) .GT. 0  ) THEN
235            DRAG(J) = ( CDS+FMASK(J)*(CDL-CDS) )
236         &                         * DENVV(J,3) * GRDSIG(NLEVxy(J,myThid))
237           ELSE
238            DRAG(J) = 0.
239         ENDIF         ENDIF
240        ENDDO        ENDDO
241    c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
242  C  C
243  C     2.2 Sensible heat flux (from clim. TS over land)  C     2.2 Sensible heat flux (from clim. TS over land)
244  C  C
# Line 256  c        CALL DUMP_WRITE2D ( NGP,1,'AUX. Line 310  c        CALL DUMP_WRITE2D ( NGP,1,'AUX.
310  C  C
311  C     2.3 Evaporation  C     2.3 Evaporation
312  C  C
313        CALL SHTORH (0,NGP,TLAND,PSA,1. _d 0,QDUMMY,RDUMMY,QSAT0(1,1))        CALL SHTORH (0,NGP,TLAND,PSA,1. _d 0,QDUMMY,RDUMMY,QSAT0(1,1),myThid)
314        CALL SHTORH (0,NGP,TSEA ,PSA,1. _d 0,QDUMMY,RDUMMY,QSAT0(1,2))        CALL SHTORH (0,NGP,TSEA ,PSA,1. _d 0,QDUMMY,RDUMMY,QSAT0(1,2),myThid)
315    
316        DO J=1,NGP        DO J=1,NGP
317  Cdj        EVAP(J,1) = CHL*DENVV(J,1)*MAX(0. _d 0,SWAV(J)*QSAT0(J,1)-Q0(J))  Cdj        EVAP(J,1) = CHL*DENVV(J,1)*MAX(0. _d 0,SWAV(J)*QSAT0(J,1)-Q0(J))
318          EVAP(J,1)=CHL*DENVV(J,1)*MIN(1. _d 0,SWAV(J))*(QSAT0(J,1)-Q0(J))  c_jmc: EVAP(J,1)=CHL*DENVV(J,1)*MIN(1. _d 0,SWAV(J))*(QSAT0(J,1)-Q0(J))
319           EVAP(J,1)=CHL*DENVV(J,1)*MAX(0. _d 0,SWAV(J)*(QSAT0(J,1)-Q0(J)))
320  cdj try new scheme : assume that the net heat flux on land is zero  cdj try new scheme : assume that the net heat flux on land is zero
321  cdj                  at all time and at all points (this is equivalent  cdj                  at all time and at all points (this is equivalent
322  cdj                  to say that the land has a zero heat capacity)  cdj                  to say that the land has a zero heat capacity)
323  cdj        EVAP(J,1) = SSR(J)-SLR(J)-SHF(J,1)  cdj        EVAP(J,1) = SSR(J)-SLR(J)-SHF(J,1)
324          EVAP(J,2) = CHS*DENVV(J,2)*(QSAT0(J,2)-Q0(J))  c_jmc:  EVAP(J,2) = CHS*DENVV(J,2)*(QSAT0(J,2)-Q0(J))
325            EVAP(J,2) = CHS*DENVV(J,2)*MAX(0. _d 0,QSAT0(J,2)-Q0(J))
326        ENDDO        ENDDO
327    
328    
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22