/[MITgcm]/MITgcm/pkg/exf/exf_wind.F
ViewVC logotype

Annotation of /MITgcm/pkg/exf/exf_wind.F

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


Revision 1.2 - (hide annotations) (download)
Tue May 30 22:47:08 2006 UTC (18 years ago) by mlosch
Branch: MAIN
Changes since 1.1: +2 -1 lines
o make the wind blow (again), if wind speed is read from file

1 mlosch 1.2 c $Header: /u/gcmpack/MITgcm/pkg/exf/exf_wind.F,v 1.1 2006/05/25 18:32:56 heimbach Exp $
2 heimbach 1.1
3     #include "EXF_OPTIONS.h"
4    
5     subroutine exf_wind(mytime, myiter, mythid)
6    
7     c ==================================================================
8     c SUBROUTINE exf_wind
9     c ==================================================================
10     c
11     c o Prepare wind speed and stress fields
12     c
13     c ==================================================================
14     c SUBROUTINE exf_wind
15     c ==================================================================
16    
17     implicit none
18    
19     c == global variables ==
20    
21     #include "EEPARAMS.h"
22     #include "SIZE.h"
23     #include "PARAMS.h"
24     #include "DYNVARS.h"
25     #include "GRID.h"
26    
27     #include "exf_param.h"
28     #include "exf_fields.h"
29     #include "exf_constants.h"
30    
31     #ifdef ALLOW_AUTODIFF_TAMC
32     #include "tamc.h"
33     #endif
34    
35     c == routine arguments ==
36    
37     integer mythid
38     integer myiter
39     _RL mytime
40    
41     #ifdef ALLOW_BULKFORMULAE
42    
43     c == local variables ==
44    
45     integer bi,bj
46     integer i,j,k
47    
48     _RL ustmp
49     _RL ustar
50    
51     c == external functions ==
52    
53     integer ilnblnk
54     external ilnblnk
55    
56     #ifndef ALLOW_ATM_WIND
57     _RL TMP1
58     _RL TMP2
59     _RL TMP3
60     _RL TMP4
61     _RL TMP5
62     #endif
63    
64     c == end of interface ==
65    
66     c-- Use atmospheric state to compute surface fluxes.
67    
68     c Loop over tiles.
69     do bj = mybylo(mythid),mybyhi(mythid)
70     do bi = mybxlo(mythid),mybxhi(mythid)
71     k = 1
72     do j = 1,sny
73     do i = 1,snx
74     c
75     c-- Initialise
76     us(i,j,bi,bj) = 0. _d 0
77     cw(i,j,bi,bj) = 0. _d 0
78     sw(i,j,bi,bj) = 0. _d 0
79     sh(i,j,bi,bj) = 0. _d 0
80     c
81     #ifdef ALLOW_ATM_WIND
82     c Wind speed and direction.
83     ustmp = uwind(i,j,bi,bj)*uwind(i,j,bi,bj) +
84     & vwind(i,j,bi,bj)*vwind(i,j,bi,bj)
85     if ( ustmp .ne. 0. _d 0 ) then
86     us(i,j,bi,bj) = sqrt(ustmp)
87     cw(i,j,bi,bj) = uwind(i,j,bi,bj)/us(i,j,bi,bj)
88     sw(i,j,bi,bj) = vwind(i,j,bi,bj)/us(i,j,bi,bj)
89     else
90     us(i,j,bi,bj) = 0. _d 0
91     cw(i,j,bi,bj) = 0. _d 0
92     sw(i,j,bi,bj) = 0. _d 0
93     endif
94     #else /* ifndef ALLOW_ATM_WIND */
95     c
96     #ifdef ALLOW_ATM_TEMP
97     c
98     c The variables us, sh and rdn have to be computed from
99     c given wind stresses inverting relationship for neutral
100     c drag coeff. cdn.
101     c The inversion is based on linear and quadratic form of
102     c cdn(umps); ustar can be directly computed from stress;
103    
104     ustmp = ustress(i,j,bi,bj)*ustress(i,j,bi,bj) +
105     & vstress(i,j,bi,bj)*vstress(i,j,bi,bj)
106     if ( ustmp .ne. 0. _d 0 ) then
107     ustar = sqrt(ustmp/atmrho)
108     cw(i,j,bi,bj) = ustress(i,j,bi,bj)/sqrt(ustmp)
109     sw(i,j,bi,bj) = vstress(i,j,bi,bj)/sqrt(ustmp)
110     else
111     ustar = 0. _d 0
112     cw(i,j,bi,bj) = 0. _d 0
113     sw(i,j,bi,bj) = 0. _d 0
114     endif
115    
116     if ( ustar .eq. 0. _d 0 ) then
117     us(i,j,bi,bj) = 0. _d 0
118     else if ( ustar .lt. ustofu11 ) then
119     tmp1 = -cquadrag_2/cquadrag_1/2
120     tmp2 = sqrt(tmp1*tmp1 + ustar*ustar/cquadrag_1)
121     us(i,j,bi,bj) = sqrt(tmp1 + tmp2)
122     else
123     tmp3 = clindrag_2/clindrag_1/3
124     tmp4 = ustar*ustar/clindrag_1/2 - tmp3**3
125     tmp5 = sqrt(ustar*ustar/clindrag_1*
126     & (ustar*ustar/clindrag_1/4 - tmp3**3))
127     us(i,j,bi,bj) = (tmp4 + tmp5)**(1/3) +
128     & tmp3**2 * (tmp4 + tmp5)**(-1/3) - tmp3
129     endif
130     c
131     #endif /* ALLOW_ATM_TEMP */
132     c
133     #endif /* ifndef ALLOW_ATM_WIND */
134    
135     c-- set lower limit
136     sh(i,j,bi,bj) = max(us(i,j,bi,bj),umin)
137    
138     c-- if wspeed available, overwrite sh,
139     c-- otherwise fill wspeed array for later use
140     if ( wspeedfile .NE. ' ' ) then
141 mlosch 1.2 us(i,j,bi,bj) = wspeed(i,j,bi,bj)
142 heimbach 1.1 sh(i,j,bi,bj) = max(wspeed(i,j,bi,bj),umin)
143     else
144     wspeed(i,j,bi,bj) = sh(i,j,bi,bj)
145     endif
146    
147     enddo
148     enddo
149     enddo
150     enddo
151    
152     #endif /* ALLOW_BULKFORMULAE */
153    
154     end

  ViewVC Help
Powered by ViewVC 1.1.22