/[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.3 - (hide annotations) (download)
Sat Jun 3 21:36:15 2006 UTC (17 years, 11 months 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 mlosch 1.3 c $Header: /u/gcmpack/MITgcm/pkg/exf/exf_wind.F,v 1.2 2006/05/30 22:47:08 mlosch 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     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 mlosch 1.3 _RL tmp1
55     _RL tmp2
56     _RL tmp3
57     _RL tmp4
58     _RL tmp5
59 heimbach 1.1
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 mlosch 1.3 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 heimbach 1.1 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 mlosch 1.2 us(i,j,bi,bj) = wspeed(i,j,bi,bj)
136 heimbach 1.1 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 mlosch 1.3 return
147 heimbach 1.1 end

  ViewVC Help
Powered by ViewVC 1.1.22