/[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.5 - (hide annotations) (download)
Thu Mar 8 08:42:24 2007 UTC (17 years, 3 months ago) by mlosch
Branch: MAIN
Changes since 1.4: +31 -29 lines
put back computation of uwind/vwind from stress from now, because it
is needed in seaice. There are doubts about the reliability of this
computation.

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

  ViewVC Help
Powered by ViewVC 1.1.22