/[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.2 by mlosch, Tue May 30 22:47:08 2006 UTC revision 1.9 by jmc, Mon May 14 19:35:39 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    
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    
 #include "EEPARAMS.h"  
22  #include "SIZE.h"  #include "SIZE.h"
23    #include "EEPARAMS.h"
24  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
 #include "GRID.h"  
25    
26  #include "exf_param.h"  #include "EXF_PARAM.h"
27  #include "exf_fields.h"  #include "EXF_FIELDS.h"
28  #include "exf_constants.h"  #include "EXF_CONSTANTS.h"
29    
30  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
31  #include "tamc.h"  #include "tamc.h"
32  #endif  #endif
33    
34  c     == routine arguments ==  C     == routine arguments ==
   
       integer mythid  
       integer myiter  
       _RL     mytime  
   
 #ifdef ALLOW_BULKFORMULAE  
35    
36  c     == local variables ==        _RL     myTime
37          INTEGER myIter
38          INTEGER myThid
39    
40        integer bi,bj  C     == external functions ==
       integer i,j,k  
41    
42        _RL     ustmp  C     == local variables ==
       _RL     ustar  
43    
44  c     == external functions ==        INTEGER bi,bj
45          INTEGER i,j
46        integer  ilnblnk        _RL     wsLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47        external ilnblnk  #ifdef ALLOW_ATM_WIND
48          _RL     wsSq
49  #ifndef ALLOW_ATM_WIND  #else
50        _RL       TMP1        _RL     usSq, recip_sqrtRhoA, ustar
51        _RL       TMP2        _RL     tmp1, tmp2, tmp3, tmp4
52        _RL       TMP3        _RL     oneThirdRL
53        _RL       TMP4        PARAMETER ( oneThirdRL = 1.d0 / 3.d0 )
       _RL       TMP5  
54  #endif  #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          k = 1  
64          do j = 1,sny  C--   Initialise
65           do i = 1,snx          DO j = 1,sNy
66  c           DO i = 1,sNx
67  c--   Initialise            wsLoc(i,j) = 0. _d 0
           us(i,j,bi,bj) = 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  #else  /* ifndef ALLOW_ATM_WIND */           ENDDO
92  c          ENDDO
93  #ifdef ALLOW_ATM_TEMP          IF ( wspeedfile .EQ. ' ' ) THEN
94  c  C-    wind-speed is not loaded from file: save local array into common block
95  c             The variables us, sh and rdn have to be computed from            DO j = 1,sNy
96  c             given wind stresses inverting relationship for neutral             DO i = 1,sNx
97  c             drag coeff. cdn.               wspeed(i,j,bi,bj) = wsLoc(i,j)
98  c             The inversion is based on linear and quadratic form of             ENDDO
99  c             cdn(umps); ustar can be directly computed from stress;            ENDDO
100            ENDIF
           ustmp = ustress(i,j,bi,bj)*ustress(i,j,bi,bj) +  
      &         vstress(i,j,bi,bj)*vstress(i,j,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  
 c          
 #endif /* ALLOW_ATM_TEMP */  
 c  
 #endif /* ifndef ALLOW_ATM_WIND */  
101    
102  c--   set lower limit  #else  /* ifndef ALLOW_ATM_WIND */
           sh(i,j,bi,bj)    = max(us(i,j,bi,bj),umin)  
103    
104  c--   if wspeed available, overwrite sh,  C--   Wind stress and direction.
105  c--   otherwise fill wspeed array for later use          DO j = 1,sNy
106            if ( wspeedfile .NE. ' ' ) then           DO i = 1,sNx
107               us(i,j,bi,bj) = wspeed(i,j,bi,bj)             IF ( stressIsOnCgrid ) THEN
108               sh(i,j,bi,bj) = max(wspeed(i,j,bi,bj),umin)               usSq = ( ustress(i,  j,bi,bj)*ustress(i  ,j,bi,bj)
109            else       &               +ustress(i+1,j,bi,bj)*ustress(i+1,j,bi,bj)
110               wspeed(i,j,bi,bj) = sh(i,j,bi,bj)       &               +vstress(i,j,  bi,bj)*vstress(i,j  ,bi,bj)
111            endif       &               +vstress(i,j+1,bi,bj)*vstress(i,j+1,bi,bj)
112         &              )*0.5 _d 0
113           enddo             ELSE
114          enddo               usSq = ustress(i,j,bi,bj)*ustress(i,j,bi,bj)
115         enddo       &             +vstress(i,j,bi,bj)*vstress(i,j,bi,bj)
116        enddo             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  #endif /* ALLOW_BULKFORMULAE */  C--   set wind-speed lower limit
180            DO j = 1,sNy
181             DO i = 1,sNx
182               sh(i,j,bi,bj) = MAX(wspeed(i,j,bi,bj),uMin)
183             ENDDO
184            ENDDO
185    
186    C--   end bi,bj loops
187           ENDDO
188          ENDDO
189    
190        end        RETURN
191          END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22