Changeset 2674 in lmdz_wrf


Ignore:
Timestamp:
Jul 12, 2019, 9:47:35 PM (5 years ago)
Author:
lfita
Message:

Adding:

  • `frontogenesis': calculation of the vector diagnostics of the frontogenesis
Location:
trunk/tools
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/diag_tools.py

    r2655 r2674  
    4242# Forcompute_fog_FRAML50: fog and visibility following Gultepe and Milbrandt, (2010)
    4343# Forcompute_front_R04: Function to compute front following Rodrigues et al.(2004)
     44# Forcompute_frontogenesis: Function to compute frontogenesis
    4445# Forcompute_potevap_orPM: Function to compute potential evapotranspiration following
    4546#   Penman-Monteith formulation implemented in ORCHIDEE
    … …  
    15761577
    15771578    return front, dt1tas, dd1wss, dt2tas, frontdims, frontvdims
     1579
     1580def Forcompute_frontogenesis(theta, ua, va, wa, press, dxs, dys, dzs, dts, dimns,    \
     1581  dimvns):
     1582    """ Function to compute frontogenesis
     1583    Forcompute_frontogenesis(ta, ua, va, wa, press, dxs, dys, dzs, dts, dimns, dimvns)
     1584      [theta]= potential air-temperature values (assuming [[t],z,y,x]) [K]
     1585      [ua]= latitudinal wind (assuming [[t],z,y,x]) [ms-1]
     1586      [va]= meridional wind (assuming [[t],z,y,x]) [ms-1]
     1587      [wa]= vertical wind (assuming [[t],z,y,x]) [ms-1]
     1588      [dxs]= grid spacing between grid points along x-axis (assuming [y,x]) [m]
     1589      [dys]= grid spacing between grid points along y-axis (assuming [y,x]) [m]
     1590      [dzs]= grid spacing between grid points along z-axis (assuming [z]) [m]
     1591      [dts]= time-step [s]
     1592      [dimns]= list of the name of the dimensions of [pa]
     1593      [dimvns]= list of the name of the variables with the values of the
     1594        dimensions of [pa]
     1595    """
     1596    fname = 'Forcompute_frontogenesis'
     1597
     1598    frontdims = dimns[:]
     1599    frontvdims = dimvns[:]
     1600
     1601    if len(theta.shape) == 4:
     1602        dx = tas.shape[3]
     1603        dy = tas.shape[2]
     1604        dz = tas.shape[1]
     1605        dt = tas.shape[0]
     1606
     1607        xdiab = np.zeros((dt,dz,dy,dx), dtype=np.float)
     1608        ydiab = np.zeros((dt,dz,dy,dx), dtype=np.float)
     1609        zdiab = np.zeros((dt,dz,dy,dx), dtype=np.float)
     1610        xdef = np.zeros((dt,dz,dy,dx), dtype=np.float)
     1611        ydef = np.zeros((dt,dz,dy,dx), dtype=np.float)
     1612        zdef = np.zeros((dt,dz,dy,dx), dtype=np.float)
     1613        xtilt = np.zeros((dt,dz,dy,dx), dtype=np.float)
     1614        ytilt = np.zeros((dt,dz,dy,dx), dtype=np.float)
     1615        zdiv = np.zeros((dt,dz,dy,dx), dtype=np.float)
     1616        f = np.zeros((3,dt,dz,dy,dx), dtype=np.float)
     1617
     1618        thetat = theta.transpose()
     1619        uat = ua.transpose()
     1620        vat = va.transpose()
     1621        wat = wa.transpose()
     1622        presst = press.transpose()
     1623        dxst = dxs.transpose()
     1624        dyst = dys.transpose()
     1625        dzst = dzs.transpose()
     1626 
     1627        pxdiab, pydiab, pzdiab, pxdef, pydef, pzdef, pxtilt, pytilt, pzdiv, pf =     \
     1628          fdin.module_fordiagnostics.compute_frontogenesis(theta=thetat, ua=uat,     \
     1629          va=vat, wa=wat, press=presst, dsx=dxst, dsy=dyst, dsz=dzst, dst=dts, d1=dx,\
     1630          d2=dy, d3=dz, d4=dt)
     1631        xdiab = pxdiab.transpose()
     1632        ydiab = pydiab.transpose()
     1633        zdiab = pzdiab.transpose()
     1634        xdef = pxdef.transpose()
     1635        ydef = pydef.transpose()
     1636        zdef = pzdef.transpose()
     1637        xtilt = pxtilt.transpose()
     1638        ytilt = pytilt.transpose()
     1639        zdiv = pzdiv.transpose()
     1640        f = pf.transpose()
     1641    else:
     1642        print errormsg
     1643        print '  ' + fname + ': rank', len(theta.shape), 'not ready !!'
     1644        print '  it only computes 4D [t,z,y,x] rank values'
     1645        quit(-1)
     1646
     1647    return xdiab, ydiab, zdiab, xdef, ydef, zdef, xtilt, ytilt, zdiv, f
    15781648
    15791649####### ###### ##### #### ### ## # END Fortran diagnostics

  • trunk/tools/diagnostics.py

    r2661 r2674  
    9090    #######
    9191availdiags = ['ACRAINTOT', 'accum', 'clt', 'cllmh', 'convini', 'deaccum', 'fog_K84', \
    92   'fog_RUC', 'front_R04', 'rhs_tas_tds', 'LMDZrh', 'mslp', 'OMEGAw', 'RAINTOT',      \
    93   'range_faces',                                                                     \
     92  'fog_RUC', 'front_R04', 'frontogenesis', 'rhs_tas_tds', 'LMDZrh', 'mslp', 'OMEGAw',\
     93  'RAINTOT', 'range_faces',                                                          \
    9494  'rvors', 'td', 'timemax', 'timeoverthres', 'turbulence', 'tws', 'uavaFROMwswd',    \
    9595  'WRFbnds', 'WRFcape_afwa', 'WRFclivi', 'WRFclwvi', 'WRF_denszint', 'WRFgeop',      \
    … …  
    169169WRFdx_compute = False
    170170WRFdx_compute = False
     171WRFdz_compute = False
    171172WRFdxdy_compute = False
    172173WRFdxdywps_compute = False
    … …  
    198199    if dnv == 'WRFdx':WRFdx_compute = True
    199200    if dnv == 'WRFdy':WRFdy_compute = True
     201    if dnv == 'WRFdz':WRFdz_compute = True
    200202    if dnv == 'WRFdxdy':WRFdxdy_compute = True
    201203    if dnv == 'WRFdxdywps':WRFdxdywps_compute = True
    … …  
    231233        if gen.searchInlist(depvars, 'WRFdx'): WRFdx_compute = True
    232234        if gen.searchInlist(depvars, 'WRFdy'): WRFdy_compute = True
     235        if gen.searchInlist(depvars, 'WRFdz'): WRFdz_compute = True
    233236        if gen.searchInlist(depvars, 'WRFdxdy'): WRFdxdy_compute = True
    234237        if gen.searchInlist(depvars, 'WRFdxdywps'): WRFdxdywps_compute = True
    … …  
    478481    Vvals = gen.variables_values('WRFds')
    479482    dictcompvars['WRFds'] = {'name': Vvals[0], 'standard_name': Vvals[1],            \
     483      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
     484
     485if WRFdz_compute:
     486    print '  ' + main + ': Retrieving dz: real distance between grid points ' +      \
     487      'from WRF as dz=PHB(k+1)+PH(k+1)-(PHB(k)+PH(k))'
     488    dimv = ncobj.variables['PH'].shape
     489
     490    dimx = dimv[3]
     491    dimy = dimv[2]
     492    dimz = dimv[1]
     493    dimt = dimv[0]
     494
     495    WRFdz = np.zeros((dimt,dimz,dimy,dimx), dtype=np.float)
     496
     497    WRFdz[:,0:dimz-1,:,:]=PH[:,1:dimz,:,:]+PHB[:,1:dimz,:,:]-(PH[:,0:dimz-1,:,:]+    \
     498      PHB[:,0:dimz-1,:,:])
     499
     500    # Attributes of the variable
     501    Vvals = gen.variables_values('WRFdz')
     502    dictcompvars['WRFdz'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
    480503      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
    481504
    … …  
    850873        ncobj.sync()
    851874
     875# frontogenesis (theta, ua, va, wa, press, dxs, dys, dzs, time)
     876    elif diagn == 'frontogenesis':
     877           
     878        var0 = ncobj.variables[depvars[0]][:]
     879        var1 = ncobj.variables[depvars[1]][:]
     880        var2 = ncobj.variables[depvars[2]][:]
     881        if depvars[3] == 'WRFpress': var3 = WRFpress
     882        else: var3 = ncobj.variables[depvars[3]]
     883        if depvars[4] == 'WRFdx': var4 = WRFdx
     884        else: var4 = ncobj.variables[depvars[4]]
     885        if depvars[5] == 'WRFdy': var5 = WRFdy
     886        else: var5 = ncobj.variables[depvars[5]]
     887        if depvars[6] == 'WRFdz': var6 = WRFdz[0,:,0,0]
     888        else: var6 = ncobj.variables[depvars[6]]
     889        if depvars[7] == 'WRFtime': var7 = WRFtime
     890        else: var7 = ncobj.variables[depvars[7]]
     891
     892        diag1, diag2, diag3, diag4, diag5, diag6, diag7, diag8, diag9, diag10,       \
     893          diagoutd, diagoutvd = diag.Forcompute_frontogenesis(var0, var1, var2, var3,\
     894          var4, var5, var6, var7, dnames,dvnames)
     895
     896        ncvar.insert_variable(ncobj, 'diabh', diag1, diagoutd, diagoutvd, newnc,     \
     897          gen.fillValueF)
     898        ovar = newnc.variables['diabh']
     899        stdn = ovar.getncattr('standard_name')
     900        ovar.setncattr('standard_name', 'x'+stdn)
     901        lngn = ovar.getncattr('long_name')
     902        ovar.setncattr('long_name', 'x component of theta '+stdn)
     903        uts = ovar.getncattr('units')
     904        ovar.setncattr('units', 'Km-1s-1')       
     905        newnc.sync()
     906        ncvar.insert_variable(ncobj, 'diabh', diag2, diagoutd, diagoutvd, newnc,     \
     907          gen.fillValueF)
     908        ovar = newnc.variables['diabh']
     909        stdn = ovar.getncattr('standard_name')
     910        ovar.setncattr('standard_name', 'y'+stdn)
     911        lngn = ovar.getncattr('long_name')
     912        ovar.setncattr('long_name', 'y component of theta '+stdn)
     913        uts = ovar.getncattr('units')
     914        ovar.setncattr('units', 'Km-1s-1')       
     915        newnc.sync()
     916        ncvar.insert_variable(ncobj, 'diabh', diag3, diagoutd, diagoutvd, newnc,     \
     917          gen.fillValueF)
     918        ovar = newnc.variables['diabh']
     919        stdn = ovar.getncattr('standard_name')
     920        ovar.setncattr('standard_name', 'z'+stdn)
     921        lngn = ovar.getncattr('long_name')
     922        ovar.setncattr('long_name', 'z component of theta '+stdn)
     923        uts = ovar.getncattr('units')
     924        ovar.setncattr('units', 'Km-1s-1')       
     925        newnc.sync()
     926
     927        ncvar.insert_variable(ncobj, 'def', diag4, diagoutd, diagoutvd, newnc,       \
     928          gen.fillValueF)
     929        ovar = newnc.variables['def']
     930        stdn = ovar.getncattr('standard_name')
     931        ovar.setncattr('standard_name', 'thetax'+stdn)
     932        lngn = ovar.getncattr('long_name')
     933        ovar.setncattr('long_name', 'x component of theta '+stdn)
     934        uts = ovar.getncattr('units')
     935        ovar.setncattr('units', 'Km-1s-1')       
     936        newnc.sync()
     937        ncvar.insert_variable(ncobj, 'def', diag5, diagoutd, diagoutvd, newnc,       \
     938          gen.fillValueF)
     939        ovar = newnc.variables['def']
     940        stdn = ovar.getncattr('standard_name')
     941        ovar.setncattr('standard_name', 'thetay'+stdn)
     942        lngn = ovar.getncattr('long_name')
     943        ovar.setncattr('long_name', 'y component of theta '+stdn)
     944        uts = ovar.getncattr('units')
     945        ovar.setncattr('units', 'Km-1s-1')       
     946        newnc.sync()
     947        ncvar.insert_variable(ncobj, 'def', diag6, diagoutd, diagoutvd, newnc,       \
     948          gen.fillValueF)
     949        ovar = newnc.variables['def']
     950        stdn = ovar.getncattr('standard_name')
     951        ovar.setncattr('standard_name', 'thetaz'+stdn)
     952        lngn = ovar.getncattr('long_name')
     953        ovar.setncattr('long_name', 'z component of theta '+stdn)
     954        uts = ovar.getncattr('units')
     955        ovar.setncattr('units', 'Km-1s-1')       
     956        newnc.sync()
     957
     958        ncvar.insert_variable(ncobj, 'tilt', diag7, diagoutd, diagoutvd, newnc,      \
     959          gen.fillValueF)
     960        ovar = newnc.variables['tilt']
     961        stdn = ovar.getncattr('standard_name')
     962        ovar.setncattr('standard_name', 'thetax'+stdn)
     963        lngn = ovar.getncattr('long_name')
     964        ovar.setncattr('long_name', 'x component of theta '+stdn)
     965        uts = ovar.getncattr('units')
     966        ovar.setncattr('units', 'Km-1s-1')       
     967        newnc.sync()
     968
     969        ncvar.insert_variable(ncobj, 'tilt', diag8, diagoutd, diagoutvd, newnc,      \
     970          gen.fillValueF)
     971        ovar = newnc.variables['tilt']
     972        stdn = ovar.getncattr('standard_name')
     973        ovar.setncattr('standard_name', 'thetay'+stdn)
     974        lngn = ovar.getncattr('long_name')
     975        ovar.setncattr('long_name', 'y component of theta '+stdn)
     976        uts = ovar.getncattr('units')
     977        ovar.setncattr('units', 'Km-1s-1')       
     978        newnc.sync()
     979
     980        ncvar.insert_variable(ncobj, 'div', diag9, diagoutd, diagoutvd, newnc,       \
     981          gen.fillValueF)
     982        ovar = newnc.variables['div']
     983        stdn = ovar.getncattr('standard_name')
     984        ovar.setncattr('standard_name', 'thetaz'+stdn)
     985        lngn = ovar.getncattr('long_name')
     986        ovar.setncattr('long_name', 'x component of theta '+stdn)
     987        uts = ovar.getncattr('units')
     988        ovar.setncattr('units', 'Km-1s-1')       
     989        newnc.sync()
     990
     991        newdim = newnc.createDimension('e', 3)
     992        newvar = newnc.createVariable('e', 'i', ('e'))
     993        newvar[:] = range(3)
     994        ncvar.newvar.basicvardef(newvar, 'component', 'axis component: x, y, z', '-')
     995        newnc.sync()
     996
     997        diagoutd = ['e'] + diagoutd
     998        diagoutvd = ['e'] + diagoutvd
     999        ncvar.insert_variable(ncobj, 'frontogenesis', diag10, diagoutd, diagoutvd,   \
     1000          newnc, gen.fillValueF)
     1001        newnc.sync()
     1002
    8521003# LMDZrh (pres, t, r)
    8531004    elif diagn == 'LMDZrh':

  • trunk/tools/module_ForDiagnostics.f90

    r2655 r2674  
    3333! compute_fog_FRAML50: fog and visibility following Gultepe and Milbrandt, (2010)
    3434! compute_front_R04d3: Subroutine to compute presence of a front following Rodrigues et al.(2004)
     35! compute_frontogenesis: Subroutine to compute the frontogenesis
    3536! compute_massvertint1D: Subroutine to vertically integrate a 1D variable in eta vertical coordinates
    3637! compute_psl_ecmwf: Compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]
    … …  
    15701571! uas: 10-m eastward wind [ms-1]
    15711572! vas: 10-m northward wind [ms-1]
    1572 ! dxs: grid spacing betweeen grid points along x-axis [m]
    1573 ! dys: grid spacing betweeen grid points along y-axis [m]
     1573! dsx: grid spacing betweeen grid points along x-axis [m]
     1574! dsy: grid spacing betweeen grid points along y-axis [m]
    15741575! ddtas: sensitivity to the thermal temporal increment [K]
    15751576! ddwss: sensitivity to the wind gradient [ms-1m-1]
    … …  
    15871588  END SUBROUTINE compute_front_R04d3
    15881589
     1590  SUBROUTINE compute_frontogenesis(d1, d2, d3, d4, theta, ua, va, wa, press, dsx, dsy, dsz, dst,      &
     1591    xdiab, ydiab, zdiab, xdef, ydef, zdef, xtilt, ytilt, zdiv, f)
     1592! Subroutine to compute the frontogenesis
     1593
     1594    IMPLICIT NONE
     1595
     1596    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
     1597    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: theta, ua, va, wa, press
     1598    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: dsx, dsy, dsz, dst
     1599    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(out)       :: xdiab, ydiab, zdiab, xdef, ydef, zdef,    &
     1600      xtilt, ytilt, zdiv, f
     1601    REAL(r_k), DIMENSION(d1,d2,d3,d4,3), INTENT(out)     :: f
     1602 
     1603! Local
     1604    INTEGER                                              :: it
     1605
     1606!!!!!!! Variables
     1607! theta: potential temperature [K]
     1608! ua: eastward wind [ms-1]
     1609! va: eastward wind [ms-1]
     1610! wa: eastward wind [ms-1]
     1611! press: pressure [Pa]
     1612! dsx: grid spacing betweeen grid points along x-axis [m]
     1613! dsy: grid spacing betweeen grid points along y-axis [m]
     1614! dsz: grid spacing betweeen grid points along z-axis [m]
     1615! dst: time-step [s]
     1616! xdiab, ydiab, zdiab: x/y/z diabatic term [Ks-1m-1]
     1617! xdef, ydef, zdef: x/y/z deformation term [Ks-1m-1]
     1618! xtilt, ytilt: x/y tilting term [Ks-1m-1]
     1619! zdiv: vertical divergence term [Ks-1m-1]
     1620! f: frontogenetical function as vector
     1621
     1622    fname = 'compute_frontogenesis'
     1623
     1624    DO it=1, d4
     1625      CALL var_Frontogenesis(dx, dy, dz, theta(:,:,:,it), ua(:,:,:,it), va(:,:,:,it), wa(:,:,:,it),   &
     1626        press(:,:,:,it), ddx, ddy, ddz, ddt, xdiab(:,:,:,it), ydiab(:,:,:,it), zdiab(:,:,:,it),       &
     1627        xdef(:,:,:,it), ydef(:,:,:,it), zdef(:,:,:,it), xtilt(:,:,:,it), ytilt(:,:,:,it),             &
     1628        zdiv(:,:,:,it), f(:,:,:,it,:))
     1629    END DO
     1630
     1631    RETURN
     1632
     1633  END SUBROUTINE compute_front_R04d3
     1634
    15891635END MODULE module_ForDiagnostics

  • trunk/tools/module_ForDiagnosticsVars.f90

    r2670 r2674  
    3535! var_fog_FRAML50: fog and visibility following Gultepe and Milbrandt, (2010)
    3636! var_front_R04: Subroutine to compute presence of a front following Rodrigues et al.(2004)
     37! var_Frontogenesis: Subroutine to compute the terms of the equation of frontogenesis
    3738! var_potevap_orPM: potential evapotranspiration following Penman-Monteith formulation implemented in ORCHIDEE
    3839! var_psl_ecmwf: sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]
    … …  
    20472048  END SUBROUTINE var_front_R04
    20482049
    2049   SUBROUTINE var_Frontogenesis(dx, dy, dz, dt, theta, ua, va, wa, press, ddx, ddy, ddz, ddt, xdiab,   \
     2050  SUBROUTINE var_Frontogenesis(dx, dy, dz, theta, ua, va, wa, press, ddx, ddy, ddz, ddt, xdiab,       &
    20502051    ydiab, zdiab, xdef, ydef, zdef, xtilt, ytilt, zdiv, f)
    20512052  ! Subroutine to compute the terms of the equation of frontogenesis
    … …  
    20542055    IMPLICIT NONE
    20552056
    2056     INTEGER, INTENT(in)                                  :: dx, dy, dz, dt
     2057    INTEGER, INTENT(in)                                  :: dx, dy, dz
    20572058    REAL(r_k), INTENT(in)                                :: ddt
    20582059    REAL(r_k), DIMENSION(dx,dy), INTENT(in)              :: ddx,ddy
    20592060    REAL(r_k), DIMENSION(dz), INTENT(in)                 :: ddz
    2060     REAL(r_k), DIMENSION(dx,dy,dz,dt), INTENT(in)        :: theta, ua, va, wa, press
    2061     REAL(r_k), DIMENSION(dx,dy,dz,dt), INTENT(out)       :: xdiab, ydiab, zdiab
    2062     REAL(r_k), DIMENSION(dx,dy,dz,dt), INTENT(out)       :: xdef, ydef, zdef
    2063     REAL(r_k), DIMENSION(dx,dy,dz,dt), INTENT(out)       :: xtilt, ytilt, zdiv
    2064     REAL(r_k), DIMENSION(dx,dy,dz,dt), INTENT(out)       :: f
     2061    REAL(r_k), DIMENSION(dx,dy,dz,), INTENT(in)          :: theta, ua, va, wa, press
     2062    REAL(r_k), DIMENSION(dx,dy,dz,), INTENT(out)         :: xdiab, ydiab, zdiab
     2063    REAL(r_k), DIMENSION(dx,dy,dz,), INTENT(out)         :: xdef, ydef, zdef
     2064    REAL(r_k), DIMENSION(dx,dy,dz,), INTENT(out)         :: xtilt, ytilt, zdiv
     2065    REAL(r_k), DIMENSION(dx,dy,dz,3), INTENT(out)        :: f
    20652066
    20662067    ! Local
    2067     INTEGER                                              :: i,j,k,l,it
     2068    INTEGER                                              :: i,j,k,l
    20682069    REAL(r_k), DIMENSION(dx,dy,dz)                       :: modthetagrad, thetadt
    20692070    REAL(r_k), DIMENSION(dx,dy,dz,3)                     :: thetagrad, thetadtgrad, uagrad, vagrad,   &
    … …  
    20812082! xtilt, ytilt: x/y tilting term [Ks-1m-1]
    20822083! zdiv: vertical divergence term [Ks-1m-1]
    2083 ! f: frontogenetical function
     2084! f: frontogenetical function as vector
    20842085
    20852086    fname = 'var_Frontogenesis'
    20862087
    2087     CALL deltat3D(dx, dy, dz, dt, theta, ddt, 'forward', thetadt)
    2088 
    2089     ! Computing components separately by time-step
    2090     DO it=1, dt
    2091       CALL gradient3D_1o(dx, dy, dz, thetadt(:,:,:,it), ddx, ddy, ddz, thetadtgrad)
    2092       CALL gradient3D_1o(dx, dy, dz, theta(:,:,:,it), ddx, ddy, ddz, thetagrad)
    2093       CALL gradient3D_1o(dx, dy, dz, ua(:,:,:,it), ddx, ddy, ddz, uagrad)
    2094       CALL gradient3D_1o(dx, dy, dz, va(:,:,:,it), ddx, ddy, ddz, vagrad)
    2095       CALL gradient3D_1o(dx, dy, dz, wa(:,:,:,it), ddx, ddy, ddz, wagrad)
    2096       CALL gradient3D_1o(dx, dy, dz, press(:,:,:,it), ddx, ddy, ddz, pressgrad)
    2097       CALL deformation3D(dx, dy, dz, thetagrad, uagrad, vagrad, thetadef)
    2098       CALL tilting3D(dx, dy, dz, thetagrad, wagrad, thetatilt)
    2099       xdiab(:,:,:,it) = thetadtgrad(:,:,:,1)
    2100 
    2101       xdef(:,:,:,it) = thetadef(:,:,:,1)
    2102       ydef(:,:,:,it) = thetadef(:,:,:,2)
    2103       zdef(:,:,:,it) = thetadef(:,:,:,3)
    2104       xtilt(:,:,:,it) = thetatilt(:,:,:,1)
    2105       ytilt(:,:,:,it) = thetatilt(:,:,:,2)
    2106       zdiv(:,:,:,it) = thetatilt(:,:,:,3)
    2107 
    2108     END DO
    2109 
     2088    ! CALL deltat3D(dx, dy, dz, dt, theta, ddt, 'forward', thetadt)
     2089
     2090    CALL gradient3D_1o(dx, dy, dz, thetadt(:,:,:), ddx, ddy, ddz, thetadtgrad)
     2091    CALL gradient3D_1o(dx, dy, dz, theta(:,:,:), ddx, ddy, ddz, thetagrad)
     2092    CALL gradient3D_1o(dx, dy, dz, ua(:,:,:), ddx, ddy, ddz, uagrad)
     2093    CALL gradient3D_1o(dx, dy, dz, va(:,:,:), ddx, ddy, ddz, vagrad)
     2094    CALL gradient3D_1o(dx, dy, dz, wa(:,:,:), ddx, ddy, ddz, wagrad)
     2095    CALL gradient3D_1o(dx, dy, dz, press(:,:,:), ddx, ddy, ddz, pressgrad)
     2096    CALL deformation3D(dx, dy, dz, thetagrad, uagrad, vagrad, thetadef)
     2097    CALL tilting3D(dx, dy, dz, thetagrad, wagrad, thetatilt)
     2098
     2099    CALL matmodule3D(dx, dy, dz, thetagrad, modthetagrad)
     2100
     2101    xdiab(:,:,:) = zeroRK
     2102    ydiab(:,:,:) = zeroRK
     2103    zdiab(:,:,:) = zeroRK
     2104
     2105    xdef(:,:,:) = thetadef(:,:,:,1)
     2106    ydef(:,:,:) = thetadef(:,:,:,2)
     2107    zdef(:,:,:) = thetadef(:,:,:,3)
     2108    xtilt(:,:,:) = thetatilt(:,:,:,1)
     2109    ytilt(:,:,:) = thetatilt(:,:,:,2)
     2110    zdiv(:,:,:) = thetatilt(:,:,:,3)
     2111
     2112    f(:,:,:,1) = 1./modthteagrad*thetagrad(:,:,:,1)*(xdiab(:,:,:) + xdef(:,:,:) +          &
     2113      xtilt(:,:,:))
     2114    f(:,:,:,2) = thetagrad(:,:,:,2)*(ydiab(:,:,:) + ydef(:,:,:) + ytilt(:,:,:))
     2115    f(:,:,:,3) = thetagrad(:,:,:,3)*(zdiab(:,:,:) + zdef(:,:,:) + zdiv(:,:,:))
    21102116
    21112117  END SUBROUTINE var_Frontogenesis

  • trunk/tools/module_scientific.f90

    r2668 r2674  
    79237923  END FUNCTION vecmodule3D
    79247924
    7925   SUBROUTINE FUNCTION matmodule3D(d1,d2,d3,mat)
     7925  SUBROUTINE FUNCTION matmodule3D(d1, d2, d3, mat, matmodule)
    79267926  ! Subroutine to compute the module of a 3D matrix with 3 components
    79277927
Note: See TracChangeset for help on using the changeset viewer.