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

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

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


Revision 1.3 - (show annotations) (download)
Sat Jun 3 21:36:15 2006 UTC (18 years ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint58l_post, checkpoint58r_post, checkpoint58n_post, checkpoint58t_post, checkpoint58h_post, checkpoint58q_post, checkpoint58j_post, checkpoint58f_post, checkpoint58i_post, checkpoint58g_post, checkpoint58o_post, checkpoint58k_post, checkpoint58s_post, checkpoint58p_post, checkpoint58m_post
Changes since 1.2: +9 -16 lines
decouple wind/stress computation from ALLOW_BULKFORMULA: always compute
either wind from stress or stress from wind

1 c $Header: /u/gcmpack/MITgcm/pkg/exf/exf_wind.F,v 1.2 2006/05/30 22:47:08 mlosch Exp $
2
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 c == local variables ==
42
43 integer bi,bj
44 integer i,j,k
45
46 _RL ustmp
47 _RL ustar
48
49 c == external functions ==
50
51 integer ilnblnk
52 external ilnblnk
53
54 _RL tmp1
55 _RL tmp2
56 _RL tmp3
57 _RL tmp4
58 _RL tmp5
59
60 c == end of interface ==
61
62 c-- Use atmospheric state to compute surface fluxes.
63
64 c Loop over tiles.
65 do bj = mybylo(mythid),mybyhi(mythid)
66 do bi = mybxlo(mythid),mybxhi(mythid)
67 k = 1
68 do j = 1,sny
69 do i = 1,snx
70 c
71 c-- Initialise
72 us(i,j,bi,bj) = 0. _d 0
73 cw(i,j,bi,bj) = 0. _d 0
74 sw(i,j,bi,bj) = 0. _d 0
75 sh(i,j,bi,bj) = 0. _d 0
76 c
77 #ifdef ALLOW_ATM_WIND
78 c Wind speed and direction.
79 ustmp = uwind(i,j,bi,bj)*uwind(i,j,bi,bj) +
80 & vwind(i,j,bi,bj)*vwind(i,j,bi,bj)
81 if ( ustmp .ne. 0. _d 0 ) then
82 us(i,j,bi,bj) = sqrt(ustmp)
83 cw(i,j,bi,bj) = uwind(i,j,bi,bj)/us(i,j,bi,bj)
84 sw(i,j,bi,bj) = vwind(i,j,bi,bj)/us(i,j,bi,bj)
85 else
86 us(i,j,bi,bj) = 0. _d 0
87 cw(i,j,bi,bj) = 0. _d 0
88 sw(i,j,bi,bj) = 0. _d 0
89 endif
90 #else /* ifndef ALLOW_ATM_WIND */
91 c
92 c The variables us, sh and rdn have to be computed from
93 c given wind stresses inverting relationship for neutral
94 c drag coeff. cdn.
95 c The inversion is based on linear and quadratic form of
96 c cdn(umps); ustar can be directly computed from stress;
97
98 ustmp = ustress(i,j,bi,bj)*ustress(i,j,bi,bj) +
99 & vstress(i,j,bi,bj)*vstress(i,j,bi,bj)
100 if ( ustmp .ne. 0. _d 0 ) then
101 ustar = sqrt(ustmp/atmrho)
102 cw(i,j,bi,bj) = ustress(i,j,bi,bj)/sqrt(ustmp)
103 sw(i,j,bi,bj) = vstress(i,j,bi,bj)/sqrt(ustmp)
104 else
105 ustar = 0. _d 0
106 cw(i,j,bi,bj) = 0. _d 0
107 sw(i,j,bi,bj) = 0. _d 0
108 endif
109
110 if ( ustar .eq. 0. _d 0 ) then
111 us(i,j,bi,bj) = 0. _d 0
112 else if ( ustar .lt. ustofu11 ) then
113 tmp1 = -cquadrag_2/cquadrag_1/2
114 tmp2 = sqrt(tmp1*tmp1 + ustar*ustar/cquadrag_1)
115 us(i,j,bi,bj) = sqrt(tmp1 + tmp2)
116 else
117 tmp3 = clindrag_2/clindrag_1/3
118 tmp4 = ustar*ustar/clindrag_1/2 - tmp3**3
119 tmp5 = sqrt(ustar*ustar/clindrag_1*
120 & (ustar*ustar/clindrag_1/4 - tmp3**3))
121 us(i,j,bi,bj) = (tmp4 + tmp5)**(1/3) +
122 & tmp3**2 * (tmp4 + tmp5)**(-1/3) - tmp3
123 endif
124 uwind(i,j,bi,bj) = us(i,j,bi,bj)*cw(i,j,bi,bj)
125 vwind(i,j,bi,bj) = us(i,j,bi,bj)*sw(i,j,bi,bj)
126 c
127 #endif /* ifndef ALLOW_ATM_WIND */
128
129 c-- set lower limit
130 sh(i,j,bi,bj) = max(us(i,j,bi,bj),umin)
131
132 c-- if wspeed available, overwrite sh,
133 c-- otherwise fill wspeed array for later use
134 if ( wspeedfile .NE. ' ' ) then
135 us(i,j,bi,bj) = wspeed(i,j,bi,bj)
136 sh(i,j,bi,bj) = max(wspeed(i,j,bi,bj),umin)
137 else
138 wspeed(i,j,bi,bj) = sh(i,j,bi,bj)
139 endif
140
141 enddo
142 enddo
143 enddo
144 enddo
145
146 return
147 end

  ViewVC Help
Powered by ViewVC 1.1.22