/[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.7 by mlosch, Thu Mar 8 15:09:52 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        _RL     ustar
49          _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) +            ustmp = 0.5*
96       &         vstress(i,j,bi,bj)*vstress(i,j,bi,bj)       &         (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            if ( ustmp .ne. 0. _d 0 ) then            if ( ustmp .ne. 0. _d 0 ) then
101               ustar = sqrt(ustmp/atmrho)               ustar = sqrt(ustmp/atmrho)
102               cw(i,j,bi,bj) = ustress(i,j,bi,bj)/sqrt(ustmp)               cw(i,j,bi,bj) = ustress(i,j,bi,bj)/sqrt(ustmp)
# Line 106  c             cdn(umps); ustar can be di Line 106  c             cdn(umps); ustar can be di
106               cw(i,j,bi,bj)    = 0. _d 0               cw(i,j,bi,bj)    = 0. _d 0
107               sw(i,j,bi,bj)    = 0. _d 0               sw(i,j,bi,bj)    = 0. _d 0
108            endif            endif
109          
110            if ( ustar .eq. 0. _d 0 ) then            if ( ustar .eq. 0. _d 0 ) then
111               us(i,j,bi,bj) = 0. _d 0               us(i,j,bi,bj) = 0. _d 0
112            else if ( ustar .lt. ustofu11 ) then            else if ( ustar .lt. ustofu11 ) then
113               tmp1 = -cquadrag_2/cquadrag_1/2               tmp1 = -cquadrag_2/cquadrag_1/2.
114               tmp2 = sqrt(tmp1*tmp1 + ustar*ustar/cquadrag_1)               tmp2 = sqrt(tmp1*tmp1 + ustar*ustar/cquadrag_1)
115               us(i,j,bi,bj) = sqrt(tmp1 + tmp2)               us(i,j,bi,bj) = sqrt(tmp1 + tmp2)
116            else            else
117               tmp3 = clindrag_2/clindrag_1/3               tmp3 = clindrag_2/clindrag_1/3.
118               tmp4 = ustar*ustar/clindrag_1/2 - tmp3**3               tmp4 = ustar*ustar/clindrag_1/2. - tmp3**3
119               tmp5 = sqrt(ustar*ustar/clindrag_1*               tmp5 = sqrt(ustar*ustar/clindrag_1*
120       &            (ustar*ustar/clindrag_1/4 - tmp3**3))       &            (ustar*ustar/clindrag_1/4. - tmp3**3))
121               us(i,j,bi,bj)   = (tmp4 + tmp5)**(1/3) +               us(i,j,bi,bj)   = (tmp4 + tmp5)**(1./3.) +
122       &            tmp3**2 * (tmp4 + tmp5)**(-1/3) - tmp3       &            tmp3**2 * (tmp4 + tmp5)**(-1./3.) - tmp3
123            endif            endif
124            uwind(i,j,bi,bj) = us(i,j,bi,bj)*cw(i,j,bi,bj)            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)            vwind(i,j,bi,bj) = us(i,j,bi,bj)*sw(i,j,bi,bj)
 c  
126  #endif /* ifndef ALLOW_ATM_WIND */  #endif /* ifndef ALLOW_ATM_WIND */
127    
128  c--   set lower limit  c--   set lower limit
129            sh(i,j,bi,bj)    = max(us(i,j,bi,bj),umin)            sh(i,j,bi,bj)    = max(us(i,j,bi,bj),umin)
130    
131  c--   if wspeed available, overwrite sh,  c--   if wspeed available, overwrite sh,
132  c--   otherwise fill wspeed array for later use  c--   otherwise fill wspeed array for later use
133            if ( wspeedfile .NE. ' ' ) then            if ( wspeedfile .NE. ' ' ) then
134               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.7
  ViewVC Help
Powered by ViewVC 1.1.22