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

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

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

revision 1.3 by mlosch, Sat Jun 3 21:36:15 2006 UTC revision 1.4 by jmc, Wed Jan 10 21:44:38 2007 UTC
# Line 1  Line 1 
1  c $Header$  C $Header$
2    C $Name$
3    
4  #include "EXF_OPTIONS.h"  #include "EXF_OPTIONS.h"
5    
# Line 18  c     ================================== Line 19  c     ==================================
19    
20  c     == global variables ==  c     == global variables ==
21    
 #include "EEPARAMS.h"  
22  #include "SIZE.h"  #include "SIZE.h"
23    #include "EEPARAMS.h"
24  #include "PARAMS.h"  #include "PARAMS.h"
25  #include "DYNVARS.h"  c #include "DYNVARS.h"
26  #include "GRID.h"  c #include "GRID.h"
27    
28  #include "exf_param.h"  #include "exf_param.h"
29  #include "exf_fields.h"  #include "exf_fields.h"
# Line 41  c     == routine arguments == Line 42  c     == routine arguments ==
42  c     == local variables ==  c     == local variables ==
43    
44        integer bi,bj        integer bi,bj
45        integer i,j,k        integer i,j
46    
47        _RL     ustmp        _RL     ustmp
48        _RL     ustar  c     _RL     ustar
49    c     _RL       tmp1, tmp2, tmp3, tmp4, tmp5
50    
51  c     == external functions ==  c     == external functions ==
52    
53        integer  ilnblnk        integer  ilnblnk
54        external ilnblnk        external ilnblnk
55    
       _RL       tmp1  
       _RL       tmp2  
       _RL       tmp3  
       _RL       tmp4  
       _RL       tmp5  
   
56  c     == end of interface ==  c     == end of interface ==
57    
58  c--   Use atmospheric state to compute surface fluxes.  c--   Use atmospheric state to compute surface fluxes.
# Line 64  c--   Use atmospheric state to compute s Line 60  c--   Use atmospheric state to compute s
60  c     Loop over tiles.  c     Loop over tiles.
61        do bj = mybylo(mythid),mybyhi(mythid)        do bj = mybylo(mythid),mybyhi(mythid)
62         do bi = mybxlo(mythid),mybxhi(mythid)         do bi = mybxlo(mythid),mybxhi(mythid)
         k = 1  
63          do j = 1,sny          do j = 1,sny
64           do i = 1,snx           do i = 1,snx
65  c  c
# Line 94  c             given wind stresses invert Line 89  c             given wind stresses invert
89  c             drag coeff. cdn.  c             drag coeff. cdn.
90  c             The inversion is based on linear and quadratic form of  c             The inversion is based on linear and quadratic form of
91  c             cdn(umps); ustar can be directly computed from stress;  c             cdn(umps); ustar can be directly computed from stress;
92    C-   no need for wind-stress inversion since everything
93    C    (ustar, ... etc ...) is derived directly from wind-stress
94    
95            ustmp = ustress(i,j,bi,bj)*ustress(i,j,bi,bj) +  c         ustmp = ustress(i,j,bi,bj)*ustress(i,j,bi,bj) +
96       &         vstress(i,j,bi,bj)*vstress(i,j,bi,bj)  c    &         vstress(i,j,bi,bj)*vstress(i,j,bi,bj)
97            if ( ustmp .ne. 0. _d 0 ) then  c         if ( ustmp .ne. 0. _d 0 ) then
98               ustar = sqrt(ustmp/atmrho)  c            ustar = sqrt(ustmp/atmrho)
99               cw(i,j,bi,bj) = ustress(i,j,bi,bj)/sqrt(ustmp)  c            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)  c            sw(i,j,bi,bj) = vstress(i,j,bi,bj)/sqrt(ustmp)
101            else  c         else
102               ustar            = 0. _d 0  c            ustar            = 0. _d 0
103               cw(i,j,bi,bj)    = 0. _d 0  c            cw(i,j,bi,bj)    = 0. _d 0
104               sw(i,j,bi,bj)    = 0. _d 0  c            sw(i,j,bi,bj)    = 0. _d 0
105            endif  c         endif
106          c
107            if ( ustar .eq. 0. _d 0 ) then  c         if ( ustar .eq. 0. _d 0 ) then
108               us(i,j,bi,bj) = 0. _d 0  c            us(i,j,bi,bj) = 0. _d 0
109            else if ( ustar .lt. ustofu11 ) then  c         else if ( ustar .lt. ustofu11 ) then
110               tmp1 = -cquadrag_2/cquadrag_1/2  c            tmp1 = -cquadrag_2/cquadrag_1/2
111               tmp2 = sqrt(tmp1*tmp1 + ustar*ustar/cquadrag_1)  c            tmp2 = sqrt(tmp1*tmp1 + ustar*ustar/cquadrag_1)
112               us(i,j,bi,bj) = sqrt(tmp1 + tmp2)  c            us(i,j,bi,bj) = sqrt(tmp1 + tmp2)
113            else  c         else
114               tmp3 = clindrag_2/clindrag_1/3  c            tmp3 = clindrag_2/clindrag_1/3
115               tmp4 = ustar*ustar/clindrag_1/2 - tmp3**3  c            tmp4 = ustar*ustar/clindrag_1/2 - tmp3**3
116               tmp5 = sqrt(ustar*ustar/clindrag_1*  c            tmp5 = sqrt(ustar*ustar/clindrag_1*
117       &            (ustar*ustar/clindrag_1/4 - tmp3**3))  c    &            (ustar*ustar/clindrag_1/4 - tmp3**3))
118               us(i,j,bi,bj)   = (tmp4 + tmp5)**(1/3) +  c            us(i,j,bi,bj)   = (tmp4 + tmp5)**(1/3) +
119       &            tmp3**2 * (tmp4 + tmp5)**(-1/3) - tmp3  c    &            tmp3**2 * (tmp4 + tmp5)**(-1/3) - tmp3
120            endif  c         endif
           uwind(i,j,bi,bj) = us(i,j,bi,bj)*cw(i,j,bi,bj)  
           vwind(i,j,bi,bj) = us(i,j,bi,bj)*sw(i,j,bi,bj)  
 c  
121  #endif /* ifndef ALLOW_ATM_WIND */  #endif /* ifndef ALLOW_ATM_WIND */
122    
123  c--   set lower limit  c--   set lower limit
124            sh(i,j,bi,bj)    = max(us(i,j,bi,bj),umin)            sh(i,j,bi,bj)    = max(us(i,j,bi,bj),umin)
125    
126  c--   if wspeed available, overwrite sh,  c--   if wspeed available, overwrite sh,
127  c--   otherwise fill wspeed array for later use  c--   otherwise fill wspeed array for later use
128            if ( wspeedfile .NE. ' ' ) then            if ( wspeedfile .NE. ' ' ) then
129               us(i,j,bi,bj) = wspeed(i,j,bi,bj)               us(i,j,bi,bj) = wspeed(i,j,bi,bj)
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22