/[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.8 - (hide annotations) (download)
Mon Apr 16 23:27:21 2007 UTC (17 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint59a, checkpoint59
Changes since 1.7: +4 -4 lines
move EXF header files from lower_case.h to UPPER_CASE.h ;
 add missing cvs Header & Name

1 jmc 1.8 C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_wind.F,v 1.7 2007/03/08 15:09:52 mlosch 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 jmc 1.8 #include "EXF_PARAM.h"
29     #include "EXF_FIELDS.h"
30     #include "EXF_CONSTANTS.h"
31 heimbach 1.1
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.6 ustmp = 0.5*
96     & (ustress(i, j,bi,bj)*ustress(i ,j,bi,bj)
97     & +ustress(i+1,j,bi,bj)*ustress(i+1,j,bi,bj)
98     & +vstress(i,j, bi,bj)*vstress(i,j ,bi,bj)
99     & +vstress(i,j+1,bi,bj)*vstress(i,j+1,bi,bj))
100 mlosch 1.5 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 mlosch 1.6 tmp1 = -cquadrag_2/cquadrag_1/2.
114 mlosch 1.5 tmp2 = sqrt(tmp1*tmp1 + ustar*ustar/cquadrag_1)
115     us(i,j,bi,bj) = sqrt(tmp1 + tmp2)
116     else
117 mlosch 1.6 tmp3 = clindrag_2/clindrag_1/3.
118     tmp4 = ustar*ustar/clindrag_1/2. - tmp3**3
119 mlosch 1.5 tmp5 = sqrt(ustar*ustar/clindrag_1*
120 mlosch 1.6 & (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 mlosch 1.5 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 heimbach 1.1 #endif /* ifndef ALLOW_ATM_WIND */
127    
128     c-- set lower limit
129     sh(i,j,bi,bj) = max(us(i,j,bi,bj),umin)
130    
131 jmc 1.4 c-- if wspeed available, overwrite sh,
132 heimbach 1.1 c-- otherwise fill wspeed array for later use
133     if ( wspeedfile .NE. ' ' ) then
134 mlosch 1.2 us(i,j,bi,bj) = wspeed(i,j,bi,bj)
135 heimbach 1.1 sh(i,j,bi,bj) = max(wspeed(i,j,bi,bj),umin)
136     else
137     wspeed(i,j,bi,bj) = sh(i,j,bi,bj)
138     endif
139    
140     enddo
141     enddo
142     enddo
143     enddo
144    
145 mlosch 1.3 return
146 heimbach 1.1 end

  ViewVC Help
Powered by ViewVC 1.1.22