/[MITgcm]/MITgcm/pkg/cheapaml/cheapaml_coare3_flux.F
ViewVC logotype

Diff of /MITgcm/pkg/cheapaml/cheapaml_coare3_flux.F

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

revision 1.15 by jmc, Sat May 25 18:09:25 2013 UTC revision 1.16 by jmc, Mon Jun 17 13:45:14 2013 UTC
# Line 15  CBOP Line 15  CBOP
15  C     !ROUTINE: CHEAPAML_COARE3_FLUX  C     !ROUTINE: CHEAPAML_COARE3_FLUX
16  C     !INTERFACE:  C     !INTERFACE:
17        SUBROUTINE CHEAPAML_COARE3_FLUX(        SUBROUTINE CHEAPAML_COARE3_FLUX(
18       I                    i,j,bi,bj, windSq,       I                    i,j,bi,bj, iceOrNot,
19         I                    tSurf, windSq,
20       O                    hf, ef, evap, Rnl, ssqt, q100, cdq, cdu,       O                    hf, ef, evap, Rnl, ssqt, q100, cdq, cdu,
21         O                    dSensdTs, dEvapdTs, dLWdTs, dQAdTs,
22       I                    myIter, myThid )       I                    myIter, myThid )
23    
24  C     !DESCRIPTION:  C     !DESCRIPTION:
# Line 28  C     === Global variables === Line 30  C     === Global variables ===
30  #include "EEPARAMS.h"  #include "EEPARAMS.h"
31  #include "PARAMS.h"  #include "PARAMS.h"
32  #include "CHEAPAML.h"  #include "CHEAPAML.h"
 #include "DYNVARS.h"  
33    
34  C     !INPUT PARAMETERS:  C     !INPUT PARAMETERS:
35  C     windSq  :: relative wind (vs surface motion) speed square  C     i, j     :: local indices of current grid-point
36  C     myIter  :: Current iteration number in simulation  C     bi, bj   :: current tile indices
37  C     myThid  :: My Thread Id number  C     iceOrNot :: 0=open water, 1=ice cover, 2=ice+snow
38        _RL windSq  C     tSurf    :: surface temperature
39    C     windSq   :: relative wind (vs surface motion) speed square
40    C     myIter   :: Current iteration number in simulation
41    C     myThid   :: My Thread Id number
42        INTEGER i,j,bi,bj        INTEGER i,j,bi,bj
43          INTEGER iceOrNot
44          _RL tSurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
45          _RL windSq(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46        INTEGER myIter, myThid        INTEGER myIter, myThid
47  C     !OUTPUT PARAMETERS:  C     !OUTPUT PARAMETERS:
48  C     cdu :: surface drag coeff (for wind stress)  C     cdu :: surface drag coeff (for wind stress)
49        _RL hf, ef, evap, Rnl, ssqt, q100, cdq, cdu        _RL hf, ef, evap, Rnl, ssqt, q100, cdq
50          _RL cdu
51    C     derivative vs surf. temp of Sensible, Evap, LW, q100
52          _RL dSensdTs, dEvapdTs, dLWdTs, dQAdTs
53  CEOP  CEOP
54    
55  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
# Line 65  C default relative humidity Line 75  C default relative humidity
75        QaR = 0.8 _d 0        QaR = 0.8 _d 0
76    
77  C sea surface temperature without skin correction  C sea surface temperature without skin correction
78        tsw=theta(i,j,1,bi,bj)  c     tsw=theta(i,j,1,bi,bj)
79        Tas=Tair(i,j,bi,bj)        tsw = tSurf(i,j)
80          Tas = Tair(i,j,bi,bj)
81    
82  C net upward long wave  C net upward long wave
83        Rnl = 0.96 _d 0*(stefan*(tsw+celsius2K)**4) !Net longwave (up = +).        Rnl = 0.96 _d 0*(stefan*(tsw+celsius2K)**4) !Net longwave (up = +).
# Line 77  C Teten''s return s air svp es in mb Line 88  C Teten''s return s air svp es in mb
88        es = es*0.98 _d 0             !reduced for salinity Kraus 1972 p. 46        es = es*0.98 _d 0             !reduced for salinity Kraus 1972 p. 46
89  C-    convert from mb to spec. humidity  kg/kg  C-    convert from mb to spec. humidity  kg/kg
90        qs = 0.62197 _d 0*es/(p0 -0.378 _d 0*es)        qs = 0.62197 _d 0*es/(p0 -0.378 _d 0*es)
91    
92        tta = Tas+celsius2K        tta = Tas+celsius2K
93  c     ttas=tta+gamma_blk*zt  c     ttas=tta+gamma_blk*zt
94  c     ttt=tta-(CheapHgrid(i,j,bi,bj) - zt)*gamma_blk  c     ttt=tta-(CheapHgrid(i,j,bi,bj) - zt)*gamma_blk
# Line 106  C Initial guesses Line 118  C Initial guesses
118        Wg = 0.5 _d 0                 !Gustiness factor initial guess        Wg = 0.5 _d 0                 !Gustiness factor initial guess
119    
120  C Air-sea differences - includes warm layer in Dt and Dq  C Air-sea differences - includes warm layer in Dt and Dq
121  c     u = (uWind(i,j,bi,bj)-uVel(i,j,1,bi,bj))**2  c     u = (uwind(i,j,bi,bj)-uVel(i,j,1,bi,bj))**2
122  c    &  + (vWind(i,j,bi,bj)-vVel(i,j,1,bi,bj))**2  c    &  + (vwind(i,j,bi,bj)-vVel(i,j,1,bi,bj))**2
123        u = windSq        u = windSq(i,j)
124        Du= SQRT(u + Wg**2 )  !include gustiness in wind spd. difference        Du= SQRT(u + Wg**2 )  !include gustiness in wind spd. difference
125        u = SQRT(u)        u = SQRT(u)
 c     Du= (u**2 + Wg**2 )**.5  !include gustiness in wind spd. difference  
126        Dt=tsw-Tas-gamma_blk*zt  !potential temperature difference.        Dt=tsw-Tas-gamma_blk*zt  !potential temperature difference.
127        Dq=qs-q        Dq=qs-q
128    
# Line 213  C *** zoq and zot fitted to results from Line 224  C *** zoq and zot fitted to results from
224        ENDDO        ENDDO
225    
226  C compute surface fluxes and other parameters  C compute surface fluxes and other parameters
 c     tau=rhoa*usr*usr*u/Du          !stress N/m2  
227        tau=rhoa*usr*usr               !stress N/m2        tau=rhoa*usr*usr               !stress N/m2
228        hf=-cpair*rhoa*usr*tsr         !sensible W/m2        hf=-cpair*rhoa*usr*tsr         !sensible W/m2
229        ef=-lath*rhoa*usr*qsr          !latent W/m2        ef=-lath*rhoa*usr*qsr          !latent W/m2
230        evap=-rhoa*usr*qsr        evap=-rhoa*usr*qsr
231        cdq = evap/Dq        cdq = evap/Dq
232        cdu = tau/Du        cdu = tau/Du
233  c     IF (.NOT.useStressOption) THEN  
 c      ustress(i,j,bi,bj)=tau*(uWind(i,j,bi,bj)-uVel(i,j,1,bi,bj))/u  
 c      vstress(i,j,bi,bj)=tau*(vWind(i,j,bi,bj)-vVel(i,j,1,bi,bj))/u  
 c     ENDIF  
234        q100=qs+qsr*(LOG(100. _d 0/zoq)-psit(100. _d 0/L))        q100=qs+qsr*(LOG(100. _d 0/zoq)-psit(100. _d 0/L))
235    
236    C--   compute derivative of surface fluxes relatice to Tsurf
237    C-    dSensdTs = -cpair*rhoa*usr*(tsr/Dt)
238          dSensdTs = cpair*rhoa*usr
239         &                *xkar/(LOG(zt/zot10)-psit(zt/L10))
240    C-    dEvapdTs  = -rhoa*usr* d/dTs(qsr)
241    C     d/dTs(qsr)= (-xkar/(LOG(zq/zoq)-psit(zq/L)) )* d/dTs(qs)
242    C     d/dTs(qs) = 0.62197 _d 0*p0/(p0 -0.378 _d 0*es)**2 * d/dTs(es)
243    C     d/dTs(es) = (0.98)* es * 17.502 _d 0 * 240.97 _d 0 / (240.97 _d 0+tsw)**2
244    C-    this simplifies (using qs) to:
245          dEvapdTs = rhoa*usr*( xkar/(LOG(zq/zoq)-psit(zq/L)) )
246         &         * qs*p0/(p0 -0.378 _d 0*es)
247    c    &         *0.98 _d 0
248         &         * 17.502 _d 0 * 240.97 _d 0 / (240.97 _d 0+tsw)**2
249    
250          if (iceornot.EQ.0) THEN
251    c       dLWdTs = 4. _d 0*ocean_emissivity*stefan*tsr*tsr*tsr
252            dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr
253          ELSEIF (iceornot.EQ.2) THEN
254    c       dLWdTs = 4. _d 0*snow_emissivity*stefan*tsr*tsr*tsr
255            dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr
256          ELSEIF (iceornot.EQ.1) THEN
257    c       dLWdTs = 4. _d 0*ice_emissivity*stefan*tsr*tsr*tsr
258            dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr
259          ENDIF
260    
261    C--   for now, ignores derivative of q100 relative to Tsurf:
262          dQAdTs = 0.
263    
264        RETURN        RETURN
265        END        END
266    
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22