Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:40:10 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
ee830839ce Andr*0001 #include "FIZHI_OPTIONS.h"
                0002       subroutine gwdrag (myid,pz,pl,ple,dpres,pkz,uz,vz,tz,qz,phis_var,
5c6d0925b4 Andr*0003      .         dudt,dvdt,dtdt,im,jm,Lm,bi,bj,istrip,npcs,imglobal)
ee830839ce Andr*0004 C***********************************************************************
                0005 C
                0006 C  PURPOSE:
                0007 C  ========
                0008 C    Driver Routine for Gravity Wave Drag
                0009 C
                0010 C  INPUT:
                0011 C  ======
                0012 C  myid  ....... Process ID
                0013 C  pz    ....... Surface Pressure [im,jm]
5c6d0925b4 Andr*0014 C  pl    ....... 3D pressure field [im,jm,Lm]
                0015 C  ple   ....... 3d pressure at model level edges [im,jm,Lm+1]
                0016 C  dpres ....... pressure difference across level [im,jm,Lm]
                0017 C  pkz   ....... pressure**kappa [im,jm,Lm]
                0018 C  uz    ....... zonal velocity [im,jm,Lm]
                0019 C  vz    ....... meridional velocity [im,jm,Lm]
                0020 C  tz    ....... temperature [im,jm,Lm]
                0021 C  qz    ....... specific humidity [im,jm,Lm]
ee830839ce Andr*0022 C  phis_var .... topography variance
                0023 C  im    ....... number of grid points in x direction
                0024 C  jm    ....... number of grid points in y direction
5c6d0925b4 Andr*0025 C  Lm    ....... number of grid points in vertical
ee830839ce Andr*0026 C  istrip ...... 'strip' length for cache size control
                0027 C  npcs  ....... number of strips
                0028 C  imglobal .... (avg) number of longitude points around the globe
                0029 C
                0030 C  INPUT/OUTPUT:
                0031 C  ============
                0032 C  dudt  ....... Updated U-Wind   Tendency including Gravity Wave Drag
                0033 C  dvdt  ....... Updated V-Wind   Tendency including Gravity Wave Drag
                0034 C  dtdt  ....... Updated Pi*Theta Tendency including Gravity Wave Drag
                0035 C
                0036 C***********************************************************************
                0037       implicit none
                0038 
                0039 c Input Variables
                0040 c ---------------
5c6d0925b4 Andr*0041       integer myid,im,jm,Lm,bi,bj,istrip,npcs,imglobal
8b150b4af9 Andr*0042       _RL pz(im,jm)
5c6d0925b4 Andr*0043       _RL pl(im,jm,Lm)
                0044       _RL ple(im,jm,Lm+1)
                0045       _RL dpres(im,jm,Lm)
                0046       _RL pkz(im,jm,Lm)
                0047       _RL uz(im,jm,Lm)
                0048       _RL vz(im,jm,Lm)
                0049       _RL tz(im,jm,Lm)
                0050       _RL qz(im,jm,Lm)
8b150b4af9 Andr*0051       _RL phis_var(im,jm)
                0052 
5c6d0925b4 Andr*0053       _RL dudt(im,jm,Lm)
                0054       _RL dvdt(im,jm,Lm)
                0055       _RL dtdt(im,jm,Lm)
ee830839ce Andr*0056 
                0057 c Local Variables
                0058 c ---------------
5c6d0925b4 Andr*0059       _RL tv(im,jm,Lm)
                0060       _RL dragu(im,jm,Lm), dragv(im,jm,Lm)
                0061       _RL dragt(im,jm,Lm) 
8b150b4af9 Andr*0062       _RL dragx(im,jm), dragy(im,jm)
                0063       _RL sumu(im,jm)
ee830839ce Andr*0064       integer nthin(im,jm),nbase(im,jm)
                0065       integer nthini, nbasei
                0066 
8b150b4af9 Andr*0067       _RL phis_std(im,jm)
ee830839ce Andr*0068 
8b150b4af9 Andr*0069       _RL std(istrip), ps(istrip)
5c6d0925b4 Andr*0070       _RL us(istrip,Lm), vs(istrip,Lm), ts(istrip,Lm)
                0071       _RL dragus(istrip,Lm), dragvs(istrip,Lm) 
8b150b4af9 Andr*0072       _RL dragxs(istrip), dragys(istrip)
924c50866c Jean*0073       _RL plstr(istrip,Lm),plestr(istrip,Lm+1),dpresstr(istrip,Lm)
ee830839ce Andr*0074       integer nthinstr(istrip),nbasestr(istrip)
                0075 
                0076       integer n,i,j,L
8b150b4af9 Andr*0077       _RL getcon, pi
                0078       _RL grav, rgas, cp, cpinv, lstar
                0079 #ifdef ALLOW_DIAGNOSTICS
                0080       logical  diagnostics_is_on
                0081       external diagnostics_is_on
                0082       _RL tmpdiag(im,jm)
                0083 #endif
ee830839ce Andr*0084 
                0085 c Initialization
                0086 c --------------
                0087       pi    = 4.0*atan(1.0)
                0088       grav  = getcon('GRAVITY')
                0089       rgas  = getcon('RGAS')
                0090       cp    = getcon('CP')
                0091       cpinv = 1.0/cp
                0092       lstar = 2*getcon('EARTH RADIUS')*cos(pi/3.0)/imglobal
                0093 
                0094 c Compute NTHIN and NBASE
                0095 c -----------------------
                0096       do j=1,jm
                0097       do i=1,im
                0098 
5c6d0925b4 Andr*0099       do nthini = 1,Lm+1
                0100        if( pz(i,j)-ple(i,j,Lm+2-nthini).gt.25. ) then
ee830839ce Andr*0101         nthin(i,j) = nthini
                0102         goto 10
                0103        endif
                0104       enddo
                0105    10 continue
5c6d0925b4 Andr*0106       do nbasei = 1,Lm+1
                0107        if( ple(i,j,Lm+2-nbasei).lt.(0.667*pz(i,j)) ) then
ee830839ce Andr*0108         nbase(i,j) = nbasei
                0109         goto 20
                0110        endif
                0111       enddo
                0112    20 continue
5c6d0925b4 Andr*0113       if( (0.667*pz(i,j))-ple(i,j,Lm+2-nbase(i,j)) .gt. 
                0114      .           ple(i,j,Lm+3-nbase(i,j))-(0.667*pz(i,j)) ) then
ee830839ce Andr*0115       nbase(i,j) = nbase(i,j)-1
                0116       endif
                0117 
                0118       enddo
                0119       enddo
5c6d0925b4 Andr*0120 
ee830839ce Andr*0121 c Compute Topography Sub-Grid Standard Deviation
5c6d0925b4 Andr*0122 c        and constrain the Maximum Value
ee830839ce Andr*0123 c ----------------------------------------------
                0124       do j=1,jm
                0125       do i=1,im
e7070f537c Cons*0126          phis_std(i,j) = min( 400.0 _d 0, sqrt( max(0.0 _d 0,
                0127      $        phis_var(i,j)) )/grav )
ee830839ce Andr*0128       enddo
                0129       enddo
                0130 
                0131 c Compute Virtual Temperatures
                0132 c ----------------------------
5c6d0925b4 Andr*0133       do L = 1,Lm
ee830839ce Andr*0134       do j = 1,jm
                0135       do i = 1,im
                0136       tv(i,j,L) = tz(i,j,L)*pkz(i,j,L)*(1.+.609*qz(i,j,L))
                0137       enddo
                0138       enddo
                0139       enddo
                0140 
5c6d0925b4 Andr*0141       do L = 1,Lm
                0142       do j = 1,jm
                0143       do i = 1,im
                0144        dragu(i,j,L) = 0.
                0145        dragv(i,j,L) = 0.
                0146        dragt(i,j,L) = 0.
                0147       enddo
                0148       enddo
                0149       enddo
                0150 
ee830839ce Andr*0151 c Call Gravity Wave Drag Paramterization on A-Grid
                0152 c ------------------------------------------------
                0153 
                0154       do n=1,npcs
                0155 
80223129cd Andr*0156       call stripit ( phis_std,std,im*jm,im*jm,istrip,1,n )
ee830839ce Andr*0157 
80223129cd Andr*0158       call stripit ( pz,ps,im*jm,im*jm,istrip,1 ,n )
                0159       call stripit ( uz,us,im*jm,im*jm,istrip,Lm,n )
                0160       call stripit ( vz,vs,im*jm,im*jm,istrip,Lm,n )
                0161       call stripit ( tv,ts,im*jm,im*jm,istrip,Lm,n )
                0162       call stripit ( pl,plstr,im*jm,im*jm,istrip,Lm,n )
924c50866c Jean*0163       call stripit ( ple,plestr,im*jm,im*jm,istrip,Lm+1,n )
80223129cd Andr*0164       call stripit ( dpres,dpresstr,im*jm,im*jm,istrip,Lm,n )
                0165       call stripitint ( nthin,nthinstr,im*jm,im*jm,istrip,1,n )
                0166       call stripitint ( nbase,nbasestr,im*jm,im*jm,istrip,1,n )
ee830839ce Andr*0167 
                0168       call GWDD ( ps,us,vs,ts,
                0169      .            dragus,dragvs,dragxs,dragys,std,
                0170      .            plstr,plestr,dpresstr,grav,rgas,cp,
5c6d0925b4 Andr*0171      .            istrip,Lm,nthinstr,nbasestr,lstar )
ee830839ce Andr*0172 
80223129cd Andr*0173       call pastit( dragus,dragu,istrip,im*jm,im*jm,Lm,n )
                0174       call pastit( dragvs,dragv,istrip,im*jm,im*jm,Lm,n )
                0175       call pastit( dragxs,dragx,istrip,im*jm,im*jm,1 ,n )
                0176       call pastit( dragys,dragy,istrip,im*jm,im*jm,1 ,n )
ee830839ce Andr*0177 
                0178       enddo
                0179 
                0180 c Add Gravity-Wave Drag to Wind and Theta Tendencies
                0181 c -------------------------------------------------- 
5c6d0925b4 Andr*0182       do L = 1,Lm
ee830839ce Andr*0183       do j = 1,jm
                0184       do i = 1,im
e7070f537c Cons*0185          dragu(i,j,L) = sign( min(0.006 _d 0,abs(dragu(i,j,L))), dragu(i
                0186      $        ,j,L) ) 
                0187          dragv(i,j,L) = sign( min(0.006 _d 0,abs(dragv(i,j,L))), dragv(i
                0188      $        ,j,L) ) 
ee830839ce Andr*0189       dragt(i,j,L) = -( uz(i,j,L)*dragu(i,j,L)+vz(i,j,L)*dragv(i,j,L) )
                0190      .                                                         *cpinv
694aa9199c Andr*0191        dudt(i,j,L) = dudt(i,j,L) + dragu(i,j,L)
                0192        dvdt(i,j,L) = dvdt(i,j,L) + dragv(i,j,L)
                0193        dtdt(i,j,L) = dtdt(i,j,L) + dragt(i,j,L)*pz(i,j)/pkz(i,j,L)
ee830839ce Andr*0194       enddo
                0195       enddo
                0196       enddo
                0197 
                0198 c Compute Diagnostics
                0199 c -------------------
8b150b4af9 Andr*0200 #ifdef ALLOW_DIAGNOSTICS
5c6d0925b4 Andr*0201       do L = 1,Lm
8b150b4af9 Andr*0202 
                0203       if(diagnostics_is_on('GWDU    ',myid) ) then
                0204        do j=1,jm
                0205        do i=1,im
                0206         tmpdiag(i,j) = dragu(i,j,L)*86400
                0207        enddo
                0208        enddo
                0209        call diagnostics_fill(tmpdiag,'GWDU    ',L,1,3,bi,bj,myid)
ee830839ce Andr*0210       endif
                0211 
8b150b4af9 Andr*0212       if(diagnostics_is_on('GWDV    ',myid) ) then
                0213        do j=1,jm
                0214        do i=1,im
                0215         tmpdiag(i,j) = dragv(i,j,L)*86400
                0216        enddo
                0217        enddo
                0218        call diagnostics_fill(tmpdiag,'GWDV    ',L,1,3,bi,bj,myid)
                0219       endif
                0220 
                0221       if(diagnostics_is_on('GWDT    ',myid) ) then
                0222        do j=1,jm
                0223        do i=1,im
                0224         tmpdiag(i,j) = dragt(i,j,L)*86400
                0225        enddo
                0226        enddo
                0227        call diagnostics_fill(tmpdiag,'GWDT    ',L,1,3,bi,bj,myid)
                0228       endif
                0229 
                0230       enddo
                0231 
ee830839ce Andr*0232 c Gravity Wave Drag at Surface (U-Wind)
                0233 c -------------------------------------
8b150b4af9 Andr*0234       if(diagnostics_is_on('GWDUS   ',myid) ) then
                0235        call diagnostics_fill(dragx,'GWDUS   ',0,1,3,bi,bj,myid)
ee830839ce Andr*0236       endif
                0237 
                0238 c Gravity Wave Drag at Surface (V-Wind)
                0239 c -------------------------------------
8b150b4af9 Andr*0240       if(diagnostics_is_on('GWDVS   ',myid) ) then
                0241        call diagnostics_fill(dragy,'GWDVS   ',0,1,3,bi,bj,myid)
ee830839ce Andr*0242       endif
                0243 
                0244 c Gravity Wave Drag at Model Top (U-Wind)
                0245 c ---------------------------------------
8b150b4af9 Andr*0246       if(diagnostics_is_on('GWDUT   ',myid) ) then
ee830839ce Andr*0247       do j = 1,jm
                0248       do i = 1,im
                0249       sumu(i,j) = 0.0
                0250       enddo
                0251       enddo
5c6d0925b4 Andr*0252       do L = 1,Lm
ee830839ce Andr*0253       do j = 1,jm
                0254       do i = 1,im
                0255       sumu(i,j) = sumu(i,j) + dragu(i,j,L)*dpres(i,j,L)/pz(i,j)
                0256       enddo
                0257       enddo
                0258       enddo
8b150b4af9 Andr*0259        do j=1,jm
                0260        do i=1,im
                0261         tmpdiag(i,j) = dragx(i,j) + sumu(i,j)*pz(i,j)/grav*100
                0262        enddo
                0263        enddo
                0264        call diagnostics_fill(tmpdiag,'GWDUT   ',0,1,3,bi,bj,myid)
ee830839ce Andr*0265       endif
                0266 
                0267 c Gravity Wave Drag at Model Top (V-Wind)
                0268 c ---------------------------------------
8b150b4af9 Andr*0269       if(diagnostics_is_on('GWDVT   ',myid) ) then
ee830839ce Andr*0270       do j = 1,jm
                0271       do i = 1,im
                0272       sumu(i,j) = 0.0
                0273       enddo
                0274       enddo
5c6d0925b4 Andr*0275       do L = 1,Lm
ee830839ce Andr*0276       do j = 1,jm
                0277       do i = 1,im
                0278       sumu(i,j) = sumu(i,j) + dragv(i,j,L)*dpres(i,j,L)/pz(i,j)
                0279       enddo
                0280       enddo
                0281       enddo
8b150b4af9 Andr*0282        do j=1,jm
                0283        do i=1,im
                0284         tmpdiag(i,j) = dragy(i,j) + sumu(i,j)*pz(i,j)/grav*100
                0285        enddo
                0286        enddo
                0287        call diagnostics_fill(tmpdiag,'GWDVT   ',0,1,3,bi,bj,myid)
ee830839ce Andr*0288       endif
8b150b4af9 Andr*0289 #endif
ee830839ce Andr*0290 
                0291       return
                0292       end
                0293       SUBROUTINE GWDD ( ps,u,v,t,dudt,dvdt,xdrag,ydrag,
                0294      .                  std,pl,ple,dpres,
5c6d0925b4 Andr*0295      .                  grav,rgas,cp,irun,Lm,nthin,nbase,lstar )
ee830839ce Andr*0296 C***********************************************************************
                0297 C
                0298 C Description:
                0299 C  ============
                0300 C    Parameterization to introduce a Gravity Wave Drag
                0301 C    due to sub-grid scale orographic forcing
                0302 C
                0303 C Input:
                0304 C  ======
                0305 C    ps ......... Surface Pressure
                0306 C    u .......... Zonal      Wind (m/sec)
                0307 C    v .......... Meridional Wind (m/sec)
                0308 C    t .......... Virtual Temperature (deg K)
                0309 C    std ........ Standard Deviation of sub-grid Orography (m)
                0310 C    ple  ....... Model pressure Edge Values
                0311 C    pl  ........ Model pressure Values
                0312 C    dpres....... Model Delta pressure Values
                0313 C    grav ....... Gravitational constant (m/sec**2)
                0314 C    rgas ....... Gas constant
                0315 C    cp ......... Specific Heat at constant pressure
                0316 C    irun ....... Number of grid-points in horizontal dimension
5c6d0925b4 Andr*0317 C    Lm ......... Number of grid-points in vertical   dimension
ee830839ce Andr*0318 C    lstar ...... Monochromatic Wavelength/(2*pi)
                0319 C
                0320 C Output:
                0321 C  =======
                0322 C    dudt ....... Zonal Acceleration due to GW Drag (m/sec**2)
                0323 C    dvdt ....... Meridional Acceleration due to GW Drag (m/sec**2)
                0324 C    xdrag ...... Zonal Surface and Base Layer Stress (Pa)
                0325 C    ydrag ...... Meridional Surface and Base Layer Stress (Pa)
                0326 C
5c6d0925b4 Andr*0327 C NOTE: Quantities computed locally in GWDD use a
                0328 C              bottom-up counting of levels
                0329 C       The fizhi code uses a top-down so all
                0330 C       Quantities that came in through the arg list
                0331 C       must use reverse vertical indexing!!!
ee830839ce Andr*0332 C***********************************************************************
                0333 
                0334       implicit none
                0335 
                0336 c Input Variables
                0337 c ---------------
5c6d0925b4 Andr*0338       integer irun,Lm
8b150b4af9 Andr*0339       _RL ps(irun)
5c6d0925b4 Andr*0340       _RL u(irun,Lm), v(irun,Lm), t(irun,Lm)
                0341       _RL dudt(irun,Lm), dvdt(irun,Lm)
8b150b4af9 Andr*0342       _RL xdrag(irun), ydrag(irun)
                0343       _RL std(irun)
5c6d0925b4 Andr*0344       _RL ple(irun,Lm+1), pl(irun,Lm), dpres(irun,Lm)
8b150b4af9 Andr*0345       _RL grav, rgas, cp
ee830839ce Andr*0346       integer nthin(irun),nbase(irun)
8b150b4af9 Andr*0347       _RL lstar
ee830839ce Andr*0348 
                0349 c Dynamic Allocation Variables
                0350 c ----------------------------
8b150b4af9 Andr*0351       _RL ubar(irun), vbar(irun), robar(irun)
                0352       _RL speed(irun), ang(irun)
5c6d0925b4 Andr*0353       _RL bv(irun,Lm)
8b150b4af9 Andr*0354       _RL nbar(irun)
ee830839ce Andr*0355 
5c6d0925b4 Andr*0356       _RL XTENS(irun,Lm+1), YTENS(irun,Lm+1)
                0357       _RL TENSIO(irun,Lm+1)
8b150b4af9 Andr*0358       _RL DRAGSF(irun)
5c6d0925b4 Andr*0359       _RL RO(irun,Lm), DZ(irun,Lm)
ee830839ce Andr*0360 
                0361       integer icrilv(irun)
                0362 
                0363 c Local Variables
                0364 c ---------------
5c6d0925b4 Andr*0365       integer  i,L
                0366       _RL a,g,agrav,akwnmb
8b150b4af9 Andr*0367       _RL gocp,roave,roiave,frsf,gstar,vai1,vai2
                0368       _RL vaisd,velco,deluu,delvv,delve2,delz,vsqua
                0369       _RL richsn,crifro,crif2,fro2,coef
ee830839ce Andr*0370 
                0371 c Initialization
                0372 c --------------
                0373       a      = 1.0
                0374       g      = 1.0
5c6d0925b4 Andr*0375       agrav  = 1.0/grav
ee830839ce Andr*0376       akwnmb = 1.0/lstar
5c6d0925b4 Andr*0377       gocp   = grav/cp
ee830839ce Andr*0378 
5c6d0925b4 Andr*0379 c Compute Atmospheric Density (with virtual temp)
                0380 c -----------------------------------------------
                0381       do l = 1,Lm
ee830839ce Andr*0382       do i = 1,irun
5c6d0925b4 Andr*0383        ro(i,L) = pl(i,Lm+1-L)/(rgas*t(i,Lm+1-L))
ee830839ce Andr*0384       enddo
                0385       enddo
                0386 
                0387 c Compute Layer Thicknesses
                0388 c -------------------------
5c6d0925b4 Andr*0389       do l = 2,Lm
ee830839ce Andr*0390       do i = 1,irun
5c6d0925b4 Andr*0391        roiave  = ( 1./ro(i,L-1) + 1./ro(i,L) )*0.5
694aa9199c Andr*0392        dz(i,L) = agrav*roiave*( pl(i,Lm+2-L)-pl(i,Lm+1-L) )
ee830839ce Andr*0393       enddo
                0394       enddo
                0395 
                0396 
5c6d0925b4 Andr*0397 c***********************************************************************
                0398 c          Surface and Base Layer Stress                               *
                0399 c***********************************************************************
ee830839ce Andr*0400 
                0401 c Definition of Surface Wind Vector
                0402 c ---------------------------------
                0403       do  i = 1,irun
5c6d0925b4 Andr*0404        robar(i) = 0.0
ee830839ce Andr*0405        ubar(i) = 0.0
                0406        vbar(i) = 0.0
                0407       enddo
                0408 
                0409       do  i = 1,irun
                0410       do  L = 1,nbase(i)-1
694aa9199c Andr*0411        robar(i) = robar(i) + ro(i,L) * (ple(i,Lm+2-L)-ple(i,Lm+1-L))
                0412        ubar(i) = ubar(i) + u(i,Lm+1-L) * (ple(i,Lm+2-L)-ple(i,Lm+1-L))
                0413        vbar(i) = vbar(i) + v(i,Lm+1-L) * (ple(i,Lm+2-L)-ple(i,Lm+1-L))
ee830839ce Andr*0414       enddo
                0415       enddo
                0416 
                0417       do  i = 1,irun
af0ffc68df Andr*0418        robar(i) = robar(i)/(ps(i)-ple(i,Lm+1-(nbase(i)-1))) * 100.0
                0419        ubar(i) = ubar(i)/(ps(i)-ple(i,Lm+1-(nbase(i)-1)))
                0420        vbar(i) = vbar(i)/(ps(i)-ple(i,Lm+1-(nbase(i)-1)))
ee830839ce Andr*0421 
5c6d0925b4 Andr*0422        speed(i) = sqrt( ubar(i)*ubar(i) + vbar(i)*vbar(i) )
                0423        ang(i) = atan2(vbar(i),ubar(i))
ee830839ce Andr*0424       enddo
                0425 
                0426 c Brunt Vaisala Frequency
                0427 c -----------------------
                0428       do i = 1,irun
5c6d0925b4 Andr*0429        do l = 2,nbase(i)
                0430         vai1 = (t(i,Lm+1-L)-t(i,Lm+2-L))/dz(i,L)+gocp
                0431         if( vai1.LT.0.0 ) then
                0432          vai1 =  0.0
                0433         endif
                0434         vai2    = 2.0*grav/( t(i,Lm+1-L)+t(i,Lm+2-L) )
                0435         vsqua   = vai1*vai2
                0436         bv(i,L) = sqrt(vsqua)
                0437        enddo
ee830839ce Andr*0438       enddo
                0439 
                0440 c Stress at the Surface Level
                0441 c ---------------------------
                0442       do i = 1,irun
5c6d0925b4 Andr*0443        nbar(i) = 0.0
ee830839ce Andr*0444       enddo
                0445       do i = 1,irun
                0446       do l = 2,nbase(i)
694aa9199c Andr*0447        nbar(i) = nbar(i) + bv(i,L)*(pl(i,Lm+2-L)-pl(i,Lm+1-L))
ee830839ce Andr*0448       enddo
                0449       enddo
                0450 
                0451       do i = 1,irun
5c6d0925b4 Andr*0452        nbar(i) = nbar(i)/(pl(i,Lm)-pl(i,Lm+1-nbase(i)))
                0453        frsf = nbar(i)*std(i)/speed(i)
                0454 
                0455        if( speed(i).eq.0.0 .or. nbar(i).eq.0.0 ) then
                0456         tensio(i,1) = 0.0
                0457        else
                0458         gstar = g*frsf*frsf/(frsf*frsf+a*a)
                0459         tensio(i,1) = gstar*(robar(i)*speed(i)*speed(i)*speed(i))
                0460      .            / (nbar(i)*lstar)
                0461        endif
ee830839ce Andr*0462 
5c6d0925b4 Andr*0463        xtens(i,1) = tensio(i,1) * cos(ang(i))
                0464        ytens(i,1) = tensio(i,1) * sin(ang(i))
                0465        dragsf(i) = tensio(i,1)
                0466        xdrag(i) = xtens(i,1)
                0467        ydrag(i) = ytens(i,1)
ee830839ce Andr*0468       enddo
                0469 
                0470 c Check for Very thin lowest layer
                0471 c --------------------------------
                0472       do i = 1,irun
5c6d0925b4 Andr*0473        if( nthin(i).gt.1 ) then
                0474         do l = 1,nthin(i)
                0475          tensio(i,L) = tensio(i,1)
                0476          xtens(i,L) = xtens(i,1)
                0477          ytens(i,L) = ytens(i,1)
                0478         enddo
                0479        endif
ee830839ce Andr*0480       enddo
                0481 
                0482 c******************************************************
                0483 c  Compute Gravity Wave Stress from NTHIN+1 to NBASE  *
                0484 c******************************************************
                0485 
                0486       do i = 1,irun
5c6d0925b4 Andr*0487        do l = nthin(i)+1,nbase(i)
ee830839ce Andr*0488 
5c6d0925b4 Andr*0489         velco = 0.5*( (u(i,Lm+1-L)*ubar(i) + v(i,Lm+1-L)*vbar(i))
                0490      .            + (u(i,Lm+2-L)*ubar(i) + v(i,Lm+2-L)*vbar(i))  )
ee830839ce Andr*0491      .      /   speed(i)
                0492 
                0493 C Convert to Newton/m**2
5c6d0925b4 Andr*0494         roave = 0.5*(ro(i,L-1)+ro(i,L)) * 100.0     
ee830839ce Andr*0495 
5c6d0925b4 Andr*0496         if( velco.le.0.0 ) then
                0497          tensio(i,L) = tensio(i,L-1)
                0498          goto 1500
                0499         endif
ee830839ce Andr*0500                     
                0501 c Froude number squared
                0502 c ---------------------
5c6d0925b4 Andr*0503         fro2 = bv(i,L)/(akwnmb*roave*velco*velco*velco)*tensio(i,L-1)
                0504         deluu = u(i,Lm+1-L)-u(i,Lm+2-L)
                0505         delvv = v(i,Lm+1-L)-v(i,Lm+2-L)
                0506         delve2 = ( deluu*deluu + delvv*delvv )
ee830839ce Andr*0507 
                0508 c Compute Richarson Number
                0509 c ------------------------
5c6d0925b4 Andr*0510         if( delve2.ne.0.0 ) then
                0511          delz = dz(i,L)
                0512          vsqua = bv(i,L)*bv(i,L)
                0513          richsn = delz*delz*vsqua/delve2
                0514         else
                0515          richsn = 99999.0
                0516         endif
                0517 
                0518         if( richsn.le.0.25 ) then
                0519          tensio(i,L) = tensio(i,L-1)
                0520          goto 1500
                0521         endif
ee830839ce Andr*0522 
                0523 c Stress in the Base Layer changes if the local Froude number
                0524 c exceeds the Critical Froude number
                0525 c ----------------------------------
5c6d0925b4 Andr*0526         crifro = 1.0 - 0.25/richsn
                0527         crif2 = crifro*crifro
32362b8fa8 Cons*0528         if( l.eq.2 ) crif2 = min(0.7 _d 0,crif2)
ee830839ce Andr*0529 
5c6d0925b4 Andr*0530         if( fro2.gt.crif2 ) then
                0531          tensio(i,L) = crif2/fro2*tensio(i,L-1)
                0532         else
                0533          tensio(i,L) = tensio(i,L-1)
                0534         endif
ee830839ce Andr*0535 
5c6d0925b4 Andr*0536 1500    continue
                0537         xtens(i,L) = tensio(i,L)*cos(ang(i))
                0538         ytens(i,L) = tensio(i,L)*sin(ang(i))
                0539 
                0540        enddo
ee830839ce Andr*0541       enddo
                0542 
                0543 c******************************************************
                0544 c    Compute Gravity Wave Stress from Base+1 to Top   *
                0545 c******************************************************
                0546 
                0547       do i = 1,irun
5c6d0925b4 Andr*0548        icrilv(i) = 0
ee830839ce Andr*0549       enddo
                0550 
                0551       do i = 1,irun
5c6d0925b4 Andr*0552        do l = nbase(i)+1,Lm+1
ee830839ce Andr*0553 
5c6d0925b4 Andr*0554         tensio(i,L) = 0.0
ee830839ce Andr*0555 
                0556 c Check for Critical Level Absorption
                0557 c -----------------------------------
5c6d0925b4 Andr*0558         if( icrilv(i).eq.1 ) goto 130
ee830839ce Andr*0559 
                0560 c Let Remaining Stress escape out the top edge of model
                0561 c -----------------------------------------------------
5c6d0925b4 Andr*0562         if( l.eq.Lm+1 ) then
                0563          tensio(i,L) = tensio(i,L-1)
                0564          goto 130
                0565         endif
ee830839ce Andr*0566 
5c6d0925b4 Andr*0567         roave = 0.5*(ro(i,L-1)+ro(i,L)) * 100.0
                0568         vai1  = (t(i,Lm+1-L)-t(i,Lm+2-L))/dz(i,L)+gocp
ee830839ce Andr*0569  
5c6d0925b4 Andr*0570         if( vai1.lt.0.0 ) then
                0571          icrilv(i)   = 1
                0572          tensio(i,L) = 0.0
                0573          goto 130
                0574         endif
                0575 
                0576         vai2  = 2.0*grav/(t(i,Lm+1-L)+t(i,Lm+2-L))
                0577         vsqua = vai1*vai2
                0578         vaisd = sqrt(vsqua)
                0579 
                0580         velco = 0.5*( (u(i,Lm+1-L)*ubar(i) + v(i,Lm+1-L)*vbar(i))
                0581      .            + (u(i,Lm+2-L)*ubar(i) + v(i,Lm+2-L)*vbar(i))  )
ee830839ce Andr*0582      .      /   speed(i)
                0583 
5c6d0925b4 Andr*0584         if( velco.lt.0.0 ) then
                0585          icrilv(i)   = 1
                0586          tensio(i,L) = 0.0
                0587          goto 130
                0588         endif
ee830839ce Andr*0589 
                0590 c Froude number squared
                0591 c ---------------------
5c6d0925b4 Andr*0592         fro2 = vaisd/(akwnmb*roave*velco*velco*velco)*tensio(i,L-1)
                0593         deluu = u(i,Lm+1-L)-u(i,Lm+2-L)
                0594         delvv = v(i,Lm+1-L)-v(i,Lm+2-L)
                0595         delve2 = ( deluu*deluu + delvv*delvv )
ee830839ce Andr*0596 
                0597 c Compute Richarson Number
                0598 c ------------------------
5c6d0925b4 Andr*0599         if( delve2.ne.0.0 ) then
                0600          delz   = dz(i,L)
                0601          richsn = delz*delz*vsqua/delve2
                0602         else
                0603          richsn = 99999.0
                0604         endif
                0605 
                0606         if( richsn.le.0.25 ) then
                0607          tensio(i,L) = 0.0
                0608          icrilv(i)   = 1
                0609          goto 130
                0610         endif
ee830839ce Andr*0611 
                0612 c Stress in Layer changes if the local Froude number
                0613 c exceeds the Critical Froude number
                0614 c ----------------------------------
5c6d0925b4 Andr*0615         crifro = 1.0 - 0.25/richsn
                0616         crif2 = crifro*crifro
                0617 
                0618         if( fro2.ge.crif2 ) then
                0619          tensio(i,L) = crif2/fro2*tensio(i,L-1)
                0620         else
                0621          tensio(i,L) = tensio(i,L-1)
                0622         endif
                0623 
                0624   130   continue
                0625         xtens(i,L) = tensio(i,L)*cos(ang(i))
                0626         ytens(i,L) = tensio(i,L)*sin(ang(i))
                0627        enddo
ee830839ce Andr*0628       enddo
                0629  
                0630 C ******************************************************
                0631 C       MOMENTUM CHANGE FOR FREE ATMOSPHERE            *
                0632 C ******************************************************
                0633  
                0634       do i = 1,irun
5c6d0925b4 Andr*0635       do l = nthin(i)+1,Lm
694aa9199c Andr*0636        coef = -grav*ps(i)/dpres(i,Lm+1-L)
5c6d0925b4 Andr*0637        dudt(i,Lm+1-L) = coef*(xtens(i,L+1)-xtens(i,L))
                0638        dvdt(i,Lm+1-L) = coef*(ytens(i,L+1)-ytens(i,L))
ee830839ce Andr*0639       enddo
                0640       enddo
                0641  
                0642 c Momentum change near the surface
                0643 c --------------------------------
                0644       do i = 1,irun
694aa9199c Andr*0645        coef = grav*ps(i)/(ple(i,Lm+1-nthin(i))-ple(i,Lm+1))
5c6d0925b4 Andr*0646        dudt(i,Lm) = coef*(xtens(i,nthin(i)+1)-xtens(i,1))
                0647        dvdt(i,Lm) = coef*(ytens(i,nthin(i)+1)-ytens(i,1))
ee830839ce Andr*0648       enddo
                0649 
                0650 c If Lowest layer is very thin, it is strapped to next layer
                0651 c ----------------------------------------------------------
                0652       do i = 1,irun
5c6d0925b4 Andr*0653        if( nthin(i).gt.1 ) then
                0654         do l = 2,nthin(i)
                0655          dudt(i,Lm+1-L) = dudt(i,Lm)
                0656          dvdt(i,Lm+1-L) = dvdt(i,Lm)
                0657         enddo
                0658        endif
ee830839ce Andr*0659       enddo
                0660 
                0661 c Convert Units to (m/sec**2)
                0662 c --------------------------- 
5c6d0925b4 Andr*0663       do l = 1,Lm
ee830839ce Andr*0664       do i = 1,irun
5c6d0925b4 Andr*0665        dudt(i,L) = - dudt(i,L)/ps(i)*0.01
                0666        dvdt(i,L) = - dvdt(i,L)/ps(i)*0.01
ee830839ce Andr*0667       enddo
                0668       enddo
                0669 
                0670       return
                0671       end
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!