/[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.8 by jmc, Mon Apr 16 23:27:21 2007 UTC revision 1.9 by jmc, Mon May 14 19:35:39 2007 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "EXF_OPTIONS.h"  #include "EXF_OPTIONS.h"
5    
6        subroutine exf_wind(mytime, myiter, mythid)        SUBROUTINE EXF_WIND( myTime, myIter, myThid )
7    
8  c     ==================================================================  C     ==================================================================
9  c     SUBROUTINE exf_wind  C     SUBROUTINE exf_wind
10  c     ==================================================================  C     ==================================================================
11  c  C
12  c     o Prepare wind speed and stress fields  C     o Prepare wind speed and stress fields
13  c  C
14  c     ==================================================================  C     ==================================================================
15  c     SUBROUTINE exf_wind  C     SUBROUTINE exf_wind
16  c     ==================================================================  C     ==================================================================
17    
18        implicit none        IMPLICIT NONE
19    
20  c     == global variables ==  C     == global variables ==
21    
22  #include "SIZE.h"  #include "SIZE.h"
23  #include "EEPARAMS.h"  #include "EEPARAMS.h"
24  #include "PARAMS.h"  #include "PARAMS.h"
 c #include "DYNVARS.h"  
 c #include "GRID.h"  
25    
26  #include "EXF_PARAM.h"  #include "EXF_PARAM.h"
27  #include "EXF_FIELDS.h"  #include "EXF_FIELDS.h"
# Line 33  c #include "GRID.h" Line 31  c #include "GRID.h"
31  #include "tamc.h"  #include "tamc.h"
32  #endif  #endif
33    
34  c     == routine arguments ==  C     == routine arguments ==
35    
36        integer mythid        _RL     myTime
37        integer myiter        INTEGER myIter
38        _RL     mytime        INTEGER myThid
39    
40  c     == local variables ==  C     == external functions ==
41    
42        integer bi,bj  C     == local variables ==
       integer i,j  
43    
44        _RL     ustmp        INTEGER bi,bj
45        _RL     ustar        INTEGER i,j
46        _RL     tmp1, tmp2, tmp3, tmp4, tmp5        _RL     wsLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47    #ifdef ALLOW_ATM_WIND
48  c     == external functions ==        _RL     wsSq
49    #else
50        integer  ilnblnk        _RL     usSq, recip_sqrtRhoA, ustar
51        external ilnblnk        _RL     tmp1, tmp2, tmp3, tmp4
52          _RL     oneThirdRL
53          PARAMETER ( oneThirdRL = 1.d0 / 3.d0 )
54    #endif
55    
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.
59    
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)
63          do j = 1,sny  
64           do i = 1,snx  C--   Initialise
65  c          DO j = 1,sNy
66  c--   Initialise           DO i = 1,sNx
67            us(i,j,bi,bj) = 0. _d 0            wsLoc(i,j) = 0. _d 0
68            cw(i,j,bi,bj) = 0. _d 0            cw(i,j,bi,bj) = 0. _d 0
69            sw(i,j,bi,bj) = 0. _d 0            sw(i,j,bi,bj) = 0. _d 0
70            sh(i,j,bi,bj) = 0. _d 0            sh(i,j,bi,bj) = 0. _d 0
71  c            wStress(i,j,bi,bj) = 0. _d 0
72             ENDDO
73            ENDDO
74    
75  #ifdef ALLOW_ATM_WIND  #ifdef ALLOW_ATM_WIND
76  c             Wind speed and direction.  
77            ustmp = uwind(i,j,bi,bj)*uwind(i,j,bi,bj) +  C--   Wind speed and direction.
78       &         vwind(i,j,bi,bj)*vwind(i,j,bi,bj)          DO j = 1,sNy
79            if ( ustmp .ne. 0. _d 0 ) then           DO i = 1,sNx
80               us(i,j,bi,bj) = sqrt(ustmp)             wsSq = uwind(i,j,bi,bj)*uwind(i,j,bi,bj)
81               cw(i,j,bi,bj) = uwind(i,j,bi,bj)/us(i,j,bi,bj)       &         + vwind(i,j,bi,bj)*vwind(i,j,bi,bj)
82               sw(i,j,bi,bj) = vwind(i,j,bi,bj)/us(i,j,bi,bj)             IF ( wsSq .NE. 0. _d 0 ) THEN
83            else               wsLoc(i,j) = SQRT(wsSq)
84               us(i,j,bi,bj) = 0. _d 0               cw(i,j,bi,bj) = uwind(i,j,bi,bj)/wsLoc(i,j)
85                 sw(i,j,bi,bj) = vwind(i,j,bi,bj)/wsLoc(i,j)
86               ELSE
87                 wsLoc(i,j) = 0. _d 0
88               cw(i,j,bi,bj) = 0. _d 0               cw(i,j,bi,bj) = 0. _d 0
89               sw(i,j,bi,bj) = 0. _d 0               sw(i,j,bi,bj) = 0. _d 0
90            endif             ENDIF
91             ENDDO
92            ENDDO
93            IF ( wspeedfile .EQ. ' ' ) THEN
94    C-    wind-speed is not loaded from file: save local array into common block
95              DO j = 1,sNy
96               DO i = 1,sNx
97                 wspeed(i,j,bi,bj) = wsLoc(i,j)
98               ENDDO
99              ENDDO
100            ENDIF
101    
102  #else  /* ifndef ALLOW_ATM_WIND */  #else  /* ifndef ALLOW_ATM_WIND */
 c  
 c             The variables us, sh and rdn have to be computed from  
 c             given wind stresses inverting relationship for neutral  
 c             drag coeff. cdn.  
 c             The inversion is based on linear and quadratic form of  
 c             cdn(umps); ustar can be directly computed from stress;  
 C-   no need for wind-stress inversion since everything  
 C    (ustar, ... etc ...) is derived directly from wind-stress  
   
           ustmp = 0.5*  
      &         (ustress(i,  j,bi,bj)*ustress(i  ,j,bi,bj)  
      &         +ustress(i+1,j,bi,bj)*ustress(i+1,j,bi,bj)  
      &         +vstress(i,j,  bi,bj)*vstress(i,j  ,bi,bj)  
      &         +vstress(i,j+1,bi,bj)*vstress(i,j+1,bi,bj))  
           if ( ustmp .ne. 0. _d 0 ) then  
              ustar = sqrt(ustmp/atmrho)  
              cw(i,j,bi,bj) = ustress(i,j,bi,bj)/sqrt(ustmp)  
              sw(i,j,bi,bj) = vstress(i,j,bi,bj)/sqrt(ustmp)  
           else  
              ustar            = 0. _d 0  
              cw(i,j,bi,bj)    = 0. _d 0  
              sw(i,j,bi,bj)    = 0. _d 0  
           endif  
   
           if ( ustar .eq. 0. _d 0 ) then  
              us(i,j,bi,bj) = 0. _d 0  
           else if ( ustar .lt. ustofu11 ) then  
              tmp1 = -cquadrag_2/cquadrag_1/2.  
              tmp2 = sqrt(tmp1*tmp1 + ustar*ustar/cquadrag_1)  
              us(i,j,bi,bj) = sqrt(tmp1 + tmp2)  
           else  
              tmp3 = clindrag_2/clindrag_1/3.  
              tmp4 = ustar*ustar/clindrag_1/2. - tmp3**3  
              tmp5 = sqrt(ustar*ustar/clindrag_1*  
      &            (ustar*ustar/clindrag_1/4. - tmp3**3))  
              us(i,j,bi,bj)   = (tmp4 + tmp5)**(1./3.) +  
      &            tmp3**2 * (tmp4 + tmp5)**(-1./3.) - tmp3  
           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)  
 #endif /* ifndef ALLOW_ATM_WIND */  
103    
104  c--   set lower limit  C--   Wind stress and direction.
105            sh(i,j,bi,bj)    = max(us(i,j,bi,bj),umin)          DO j = 1,sNy
106             DO i = 1,sNx
107               IF ( stressIsOnCgrid ) THEN
108                 usSq = ( ustress(i,  j,bi,bj)*ustress(i  ,j,bi,bj)
109         &               +ustress(i+1,j,bi,bj)*ustress(i+1,j,bi,bj)
110         &               +vstress(i,j,  bi,bj)*vstress(i,j  ,bi,bj)
111         &               +vstress(i,j+1,bi,bj)*vstress(i,j+1,bi,bj)
112         &              )*0.5 _d 0
113               ELSE
114                 usSq = ustress(i,j,bi,bj)*ustress(i,j,bi,bj)
115         &             +vstress(i,j,bi,bj)*vstress(i,j,bi,bj)
116               ENDIF
117               IF ( usSq .NE. 0. _d 0 ) THEN
118                 wStress(i,j,bi,bj) = SQRT(usSq)
119    c            ustar = SQRT(usSq/atmrho)
120                 cw(i,j,bi,bj) = ustress(i,j,bi,bj)/wStress(i,j,bi,bj)
121                 sw(i,j,bi,bj) = vstress(i,j,bi,bj)/wStress(i,j,bi,bj)
122               ELSE
123                 wStress(i,j,bi,bj) = 0. _d 0
124                 cw(i,j,bi,bj)      = 0. _d 0
125                 sw(i,j,bi,bj)      = 0. _d 0
126               ENDIF
127             ENDDO
128            ENDDO
129    
130            IF ( wspeedfile .EQ. ' ' ) THEN
131    C--   wspeed is not loaded ; devive wind-speed by inversion of
132    C     wind-stress=fct(wind-speed) relation:
133    C             The variables us, sh and rdn have to be computed from
134    C             given wind stresses inverting relationship for neutral
135    C             drag coeff. cdn.
136    C             The inversion is based on linear and quadratic form of
137    C             cdn(umps); ustar can be directly computed from stress;
138             recip_sqrtRhoA = 1. _d 0 / SQRT(atmrho)
139             DO j = 1,sNy
140              DO i = 1,sNx
141                ustar = wStress(i,j,bi,bj)*recip_sqrtRhoA
142                IF ( ustar .EQ. 0. _d 0 ) THEN
143                 wsLoc(i,j) = 0. _d 0
144                ELSE IF ( ustar .LT. ustofu11 ) THEN
145                 tmp1 = -cquadrag_2/cquadrag_1*exf_half
146                 tmp2 = SQRT(tmp1*tmp1 + ustar*ustar/cquadrag_1)
147                 wsLoc(i,j) = SQRT(tmp1 + tmp2)
148                ELSE
149                 tmp1 = clindrag_2/clindrag_1*oneThirdRL
150                 tmp2 = ustar*ustar/clindrag_1*exf_half
151         &            - tmp1*tmp1*tmp1
152                 tmp3 = SQRT( ustar*ustar/clindrag_1*
153         &            (ustar*ustar/clindrag_1*0.25 _d 0 - tmp1*tmp1*tmp1 )
154         &                  )
155                 tmp4 = (tmp2 + tmp3)**oneThirdRL
156                 wsLoc(i,j) = tmp4 + tmp1*tmp1 / tmp4 - tmp1
157    c            wsLoc(i,j) = (tmp2 + tmp3)**oneThirdRL +
158    c    &            tmp1*tmp1 * (tmp2 + tmp3)**(-oneThirdRL) - tmp1
159                ENDIF
160              ENDDO
161             ENDDO
162    C-    save local array wind-speed to common block
163             DO j = 1,sNy
164              DO i = 1,sNx
165                wspeed(i,j,bi,bj) = wsLoc(i,j)
166              ENDDO
167             ENDDO
168            ENDIF
169    
170    C--   infer wind field from wind-speed & wind-stress direction
171            DO j = 1,sNy
172             DO i = 1,sNx
173               uwind(i,j,bi,bj) = wspeed(i,j,bi,bj)*cw(i,j,bi,bj)
174               vwind(i,j,bi,bj) = wspeed(i,j,bi,bj)*sw(i,j,bi,bj)
175             ENDDO
176            ENDDO
177    #endif /* ifndef ALLOW_ATM_WIND */
178    
179  c--   if wspeed available, overwrite sh,  C--   set wind-speed lower limit
180  c--   otherwise fill wspeed array for later use          DO j = 1,sNy
181            if ( wspeedfile .NE. ' ' ) then           DO i = 1,sNx
182               us(i,j,bi,bj) = wspeed(i,j,bi,bj)             sh(i,j,bi,bj) = MAX(wspeed(i,j,bi,bj),uMin)
183               sh(i,j,bi,bj) = max(wspeed(i,j,bi,bj),umin)           ENDDO
184            else          ENDDO
185               wspeed(i,j,bi,bj) = sh(i,j,bi,bj)  
186            endif  C--   end bi,bj loops
187           ENDDO
188           enddo        ENDDO
         enddo  
        enddo  
       enddo  
189    
190        return        RETURN
191        end        END
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22