Changeset 2274 in lmdz_wrf


Ignore:
Timestamp:
Dec 28, 2018, 3:35:45 PM (6 years ago)
Author:
lfita
Message:

Adding:

  • `WRFbnds': grid cell bounds from XLONG_U, XLAT_U, XLONG_V, XLAT_V as intersection of their related parallels and meridians
Location:
trunk/tools
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/diag_tools.py

    r2260 r2274  
    3030# compute_WRFuasvas: Fucntion to compute geographical rotated WRF 2-meter winds
    3131# derivate_centered: Function to compute the centered derivate of a given field
     32# Forcompute_cellbnds: Function to compute cellboundaries using wind-staggered lon,
     33#   lats as intersection of their related parallels and meridians
    3234# Forcompute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction
    3335#   following newmicro.F90 from LMDZ via Fortran subroutine
    … …  
    525527
    526528    return cllmh
     529
     530def Forcompute_cellbnds(ulon, ulat, vlon, vlat, dimns, dimvns):
     531    """ Function to compute cellboundaries using wind-staggered lon, lats as
     532      intersection of their related parallels and meridians
     533    compute_cellbnds(ulon, ulat, vlon, vlat, dimns, dimvns)
     534      [ulon]= x-staggered longitudes (assuming [y,x+1])
     535      [ulat]= x-staggered latitudes (assuming [y,x+1])
     536      [vlon]= y-staggered longitudes (assuming [y+1,x])
     537      [vlat]= y-staggered latitudes (assuming [y+1,x])
     538      [dimns]= list of the name of the dimensions of [cldfra]
     539      [dimvns]= list of the name of the variables with the values of the
     540        dimensions of [ulon]
     541    """
     542    fname = 'Forcompute_cellbnds'
     543
     544    dims = dimns[:]
     545    vdims = dimvns[:]
     546
     547    if len(ulon.shape) == 2:
     548        sdx = ulon.shape[1]
     549        dy = ulon.shape[0]
     550        dx = vlon.shape[1]
     551        sdy = vlon.shape[0]
     552
     553        ulont = ulon.transpose()
     554        ulatt = ulat.transpose()
     555        vlont = vlon.transpose()
     556        vlatt = vlat.transpose()
     557
     558        xbndst, ybndst = fdin.module_fordiagnostics.compute_cellbnds(ulon=ulont,     \
     559          ulat=ulatt, vlon=vlont, vlat=vlatt, dx=dx, dy=dy, sdx=sdx, sdy=sdy)
     560    else:
     561        print errormsg
     562        print '  ' + fname + ": wrong rank of variables !!"
     563        print '    2D matrices are expected and its found instead'
     564        print '    ulon shape:', ulon.shape
     565        print '    ulat shape:', ulat.shape
     566        print '    vlon shape:', vlon.shape
     567        print '    vlat shape:', vlat.shape
     568        quit(-1)
     569   
     570    xbnds = xbndst.transpose()
     571    ybnds = ybndst.transpose()
     572
     573    return xbnds, ybnds, dims, vdims
    527574
    528575def Forcompute_cllmh(cldfra, pres, dimns, dimvns):

  • trunk/tools/diagnostics.inf

    r2215 r2274  
    2222hurs, TSrhs, psfc@t@q
    2323hurs, WRFrhs, PSFC@T2@Q2
     24lat_bnds, WRFbnds, XLONG_U@XLAT_U@XLONG_V@XLAT_V
     25lon_bnds, WRFbnds, XLONG_U@XLAT_U@XLONG_V@XLAT_V
    2426mrso, WRFmrso,  SMOIS@DZS
    2527mrsos, WRFmrsos,  SMOIS@DZS

  • trunk/tools/diagnostics.py

    r2260 r2274  
    9292  'fog_RUC', 'LMDZrh', 'mslp', 'OMEGAw', 'RAINTOT', 'range_faces',                   \
    9393  'rvors', 'td', 'timemax', 'timeoverthres', 'turbulence', 'uavaFROMwswd',           \
    94   'WRFcape_afwa', 'WRFclivi', 'WRFclwvi', 'WRF_denszint', 'WRFgeop',                 \
     94  'WRFbnds', 'WRFcape_afwa', 'WRFclivi', 'WRFclwvi', 'WRF_denszint', 'WRFgeop',      \
    9595  'WRFmrso', 'WRFmrsos', 'WRFpotevap_orPM', 'WRFp', 'WRFpsl_ecmwf',                  \
    9696  'WRFpsl_ptarget', 'WRFrvors', 'WRFslw', 'ws', 'wds', 'wss', 'WRFheight',           \
    … …  
    103103NONcheckingvars = ['accum', 'cllmh', 'deaccum', 'face', 'LONLATdxdy', 'TSrhs',       \
    104104  'TStd', 'TSwds', 'TSwss',                                                          \
    105   'WRFbils',                                                                         \
     105  'WRFbils',  'WRFbnds',                                                             \
    106106  'WRFclivi', 'WRFclwvi', 'WRFdens', 'WRFdx', 'WRFdxdy', 'WRFdxdywps', 'WRFdy',      \
    107107  'WRFgeop', 'WRFp', 'WRFtd',                                                        \
    … …  
    891891        ncvar.set_attributek(ovar, 'h_valley_newrange', hvalleyrange, 'F')
    892892
     893# cell_bnds: grid cell bounds from XLONG_U, XLAT_U, XLONG_V, XLAT_V as intersection
     894#   of their related parallels and meridians
     895    elif diagn == 'WRFbnds':
     896           
     897        var0 = ncobj.variables[depvars[0]][0,:,:]
     898        var1 = ncobj.variables[depvars[1]][0,:,:]
     899        var2 = ncobj.variables[depvars[2]][0,:,:]
     900        var3 = ncobj.variables[depvars[3]][0,:,:]
     901
     902        dnamesvar = []
     903        dnamesvar.append(ncobj.variables[depvars[2]].dimensions[2])
     904        dnamesvar.append(ncobj.variables[depvars[0]].dimensions[1])
     905        dnamesvar.append('bnds')
     906        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
     907
     908        cellbndsx, cellbndsy, diagoutd, diagoutvd = diag.Forcompute_cellbnds(var0,   \
     909          var1, var2, var3, dnamesvar, dvnamesvar)
     910
     911        # Removing the nonChecking variable-dimensions from the initial list
     912        varsadd = []
     913        diagoutvd = list(dvnames)
     914        for nonvd in NONchkvardims:
     915            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
     916            varsadd.append(nonvd)
     917        # creation of bounds dimension
     918        newdim = newnc.createDimension('bnds', 4)
     919
     920        ncvar.insert_variable(ncobj, 'lon_bnds', cellbndsx, diagoutd, diagoutvd, newnc)
     921        ncvar.insert_variable(ncobj, 'lat_bnds', cellbndsy, diagoutd, diagoutvd, newnc)
     922
    893923# mrso: total soil moisture SMOIS, DZS
    894924    elif diagn == 'WRFmrso':

  • trunk/tools/module_ForDiagnostics.f90

    r2260 r2274  
    1616! compute_cape_afwa4D: Subroutine to use WRF phys/module_diag_afwa.F `buyoancy' subroutine to compute
    1717!   CAPE, CIN, ZLFC, PLFC, LI
     18! compute_cellbnds: Subroutine to compute cellboundaries using wind-staggered lon, lats as
     19!   intersection of their related parallels and meridians
    1820! compute_cllmh4D3: Computation of low, medium and high cloudiness from a 4D CLDFRA and pressure being
    1921!   3rd dimension the z-dim
    … …  
    10511053  END SUBROUTINE compute_range_faces
    10521054
     1055  SUBROUTINE compute_cellbnds(dx, dy, sdx, sdy, ulon, ulat, vlon, vlat, xbnds, ybnds)
     1056! Subroutine to compute cellboundaries using wind-staggered lon, lats as intersection of their related
     1057!   parallels and meridians
     1058
     1059    IMPLICIT NONE
     1060
     1061    INTEGER, INTENT(in)                                  :: dx, dy, sdx, sdy
     1062    REAL(r_k), DIMENSION(sdx, dy), INTENT(in)            :: ulon, ulat
     1063    REAL(r_k), DIMENSION(dx, sdy), INTENT(in)            :: vlon, vlat
     1064    REAL(r_k), DIMENSION(dx, dy, 4), INTENT(out)         :: xbnds, ybnds
     1065
     1066! Local
     1067    INTEGER                                              :: i,j,iv
     1068    INTEGER                                              :: ix,ex,iy,ey
     1069    CHARACTER(len=2), DIMENSION(4)                       :: Svertex
     1070    INTEGER, DIMENSION(4,2,2,2)                          :: indices
     1071    REAL(r_k), DIMENSION(2)                              :: ptintsct
     1072    REAL(r_k), DIMENSION(2,2)                            :: merid, paral
     1073    LOGICAL                                              :: intsct
     1074
     1075!!!!!!! Variables
     1076! dx, dy: un-staggered dimensions
     1077! sdx, sdy: staggered dimensions
     1078! ulon, ulat: x-wind staggered longitudes and latitudes
     1079! vlon, vlat: y-wind staggered longitudes and latitudes
     1080! xbnds, ybnds: x and y cell boundaries
     1081
     1082    fname = 'compute_cellbnds'
     1083
     1084    ! Indices to use indices[SW/NW/NE/SE, m/p, x/y, i/e]
     1085    Svertex = (/ 'SW', 'NW', 'NE', 'SE' /)
     1086
     1087    ! SW
     1088    indices(1,1,1,1) = 0
     1089    indices(1,1,1,2) = 0
     1090    indices(1,1,2,1) = -1
     1091    indices(1,1,2,2) = 0
     1092    indices(1,2,1,1) = -1
     1093    indices(1,2,1,2) = 0
     1094    indices(1,2,2,1) = -1
     1095    indices(1,2,2,2) = -1
     1096    ! NW
     1097    indices(2,1,1,1) = 0
     1098    indices(2,1,1,2) = 0
     1099    indices(2,1,2,1) = 0
     1100    indices(2,1,2,2) = 1
     1101    indices(2,2,1,1) = -1
     1102    indices(2,2,1,2) = 0
     1103    indices(2,2,2,1) = 1
     1104    indices(2,2,2,2) = 1
     1105    ! NE
     1106    indices(3,1,1,1) = 1
     1107    indices(3,1,1,2) = 1
     1108    indices(3,1,2,1) = 0
     1109    indices(3,1,2,2) = 1
     1110    indices(3,2,1,1) = 0
     1111    indices(3,2,1,2) = 1
     1112    indices(3,2,2,1) = 1
     1113    indices(3,2,2,2) = 1
     1114    ! SE
     1115    indices(4,1,1,1) = 1
     1116    indices(4,1,1,2) = 1
     1117    indices(4,1,2,1) = -1
     1118    indices(4,1,2,2) = 0
     1119    indices(4,2,1,1) = 0
     1120    indices(4,2,1,2) = 1
     1121    indices(4,2,2,1) = -1
     1122    indices(4,2,2,2) = -1
     1123
     1124    DO i=1,dx
     1125      DO j=1,dy
     1126        DO iv=1,4
     1127
     1128          ix = MAX(i+indices(iv,1,1,1),1)
     1129          !ex = MIN(i+indices(iv,1,1,2),dx)
     1130          ex = i+indices(iv,1,1,2)
     1131          iy = MAX(j+indices(iv,1,2,1),1)
     1132          ey = MIN(j+indices(iv,1,2,2),dy)
     1133 
     1134          merid(1,1) = ulon(ix,iy)
     1135          merid(1,2) = ulat(ix,iy)
     1136          merid(2,1) = ulon(ex,ey)
     1137          merid(2,2) = ulat(ex,ey)
     1138
     1139          ix = MAX(i+indices(iv,2,1,1),1)
     1140          ex = MIN(i+indices(iv,2,1,2),dx)
     1141          iy = MAX(j+indices(iv,2,2,1),1)
     1142          !ey = MIN(i+indices(iv,2,2,2),dy)
     1143          ey = j+indices(iv,2,2,2)
     1144          paral(1,1) = vlon(ix,iy)
     1145          paral(1,2) = vlat(ix,iy)
     1146          paral(2,1) = vlon(ex,ey)
     1147          paral(2,2) = vlat(ex,ey)
     1148
     1149          CALL intersection_2Dlines(merid, paral, intsct, ptintsct)
     1150          IF (.NOT.intsct) THEN
     1151            msg = 'not interection found for ' // Svertex(iv) // ' vertex'
     1152            CALL ErrMsg(msg, fname, -1)
     1153          END IF
     1154          xbnds(i,j,iv) = ptintsct(1)
     1155          ybnds(i,j,iv) = ptintsct(2)
     1156        END DO
     1157      END DO
     1158    END DO
     1159
     1160  END SUBROUTINE
     1161
    10531162  SUBROUTINE compute_Koeppen_Geiger_climates(dx, dy, dt, pr, tas, climatesI, climatesS, climlegend)
    10541163  ! Subroutine to compute the Koeppen-Geiger climates after:

  • trunk/tools/module_scientific.f90

    r2271 r2274  
    41404140    IF (lineA(2,1) /= lineA(1,1)) THEN
    41414141      segmentA(2) = (lineA(2,2)-lineA(1,2))/(lineA(2,1)-lineA(1,1))
    4142       IF ( (lineA(1,1)*segmentA(2) - lineA(1,2)) /= (lineA(2,1)*segmentA(2) - lineA(2,2)) ) THEN
    4143         PRINT *,'A = y1 - x2*B = ', lineA(1,2) - lineA(1,1)*segmentA(2)
    4144         PRINT *,'A = y2 - x2*B = ', lineA(2,2) - lineA(2,1)*segmentA(2)
    4145         msg = 'Wrong calculation of parameter A, for lineA'
    4146         CALL ErrMSg(msg, fname, -1)
    4147       END IF
     4142      ! This might be to ask too much... ?
     4143      !IF ( (lineA(1,1)*segmentA(2) - lineA(1,2)) /= (lineA(2,1)*segmentA(2) - lineA(2,2)) ) THEN
     4144      !  PRINT *,'A = y1 - x2*B = ', lineA(1,2) - lineA(1,1)*segmentA(2)
     4145      !  PRINT *,'A = y2 - x2*B = ', lineA(2,2) - lineA(2,1)*segmentA(2)
     4146      !  msg = 'Wrong calculation of parameter A, for lineA'
     4147      !  CALL ErrMSg(msg, fname, -1)
     4148      !END IF
    41484149      segmentA(1) = lineA(1,2) - lineA(1,1)*segmentA(2)
    41494150      a11 = segmentA(2)
    … …  
    41584159    IF (lineB(2,1) /= lineB(1,1)) THEN
    41594160      segmentB(2) = (lineB(2,2)-lineB(1,2))/(lineB(2,1)-lineB(1,1))
    4160       IF ( (lineB(1,1)*segmentB(2) - lineB(1,2)) /= (lineB(2,1)*segmentB(2) - lineB(2,2)) ) THEN
    4161         PRINT *,'A = x1*B -y1 = ', lineB(1,1)*segmentB(2) - lineB(1,2)
    4162         PRINT *,'A = x2*B -y2 = ', lineB(2,1)*segmentB(2) - lineB(2,2)
    4163         msg = 'Wrong calculation of parameter A, for lineB'
    4164         CALL ErrMSg(msg, fname, -1)
    4165       END IF
     4161      ! This might be to ask too much... ?
     4162      !IF ( (lineB(1,1)*segmentB(2) - lineB(1,2)) /= (lineB(2,1)*segmentB(2) - lineB(2,2)) ) THEN
     4163      !  PRINT *,'A = x1*B -y1 = ', lineB(1,1)*segmentB(2) - lineB(1,2)
     4164      !  PRINT *,'A = x2*B -y2 = ', lineB(2,1)*segmentB(2) - lineB(2,2)
     4165      !  msg = 'Wrong calculation of parameter A, for lineB'
     4166      !  CALL ErrMSg(msg, fname, -1)
     4167      !END IF
    41664168      segmentB(1) = lineB(1,2) - lineB(1,1)*segmentB(2)
    41674169      a21 = segmentB(2)

  • trunk/tools/variables_values.dat

    r2249 r2274  
    259259lai, lai, lai, 0., 8., Leaf|Area|Index, -, Greens, $lai$, lai
    260260LAI, lai, lai, 0., 8., Leaf|Area|Index, -, Greens, $lai$, lai
    261 lat, lat, latitude, -90., 90., latitude, degrees_North, seismic, $lat$, lat
    262 LAT, lat, latitude, -90., 90., latitude, degrees_North, seismic, $lat$, lat
    263 south_north, lat, latitude, -90., 90., latitude, degrees_North, seismic, $lat$, lat
    264 XLAT, lat, latitude, -90., 90., latitude, degrees_North, seismic, $lat$, lat
    265 XLAT_M, lat, latitude, -90., 90., latitude, degrees_North, seismic, $lat$, lat
    266 latitude, lat, latitude, -90., 90., latitude, degrees_North, seismic, $lat$, lat
    267 lat, lat, latitude, -90., 90., latitude, degrees_North, seismic, $lat$, lat
    268 nav_lat, lat, latitude, -90., 90., latitude, degrees_North, seismic, $lat$, lat
     261lat, lat, latitude, -90., 90., latitude, degrees_north, seismic, $lat$, lat
     262LAT, lat, latitude, -90., 90., latitude, degrees_north, seismic, $lat$, lat
     263south_north, lat, latitude, -90., 90., latitude, degrees_north, seismic, $lat$, lat
     264XLAT, lat, latitude, -90., 90., latitude, degrees_north, seismic, $lat$, lat
     265XLAT_M, lat, latitude, -90., 90., latitude, degrees_north, seismic, $lat$, lat
     266latitude, lat, latitude, -90., 90., latitude, degrees_north, seismic, $lat$, lat
     267lat, lat, latitude, -90., 90., latitude, degrees_north, seismic, $lat$, lat
     268nav_lat, lat, latitude, -90., 90., latitude, degrees_north, seismic, $lat$, lat
     269lat_bnds, lat_bnds, latitude_bnds, -90., 90., latitude_bounds, degrees_north, seismic, $lat_{bnds}$, lat_bnds
    269270lco, lco, land_use, 0, 24, Land|use, 1, ranibow, $lco$, lco
    270271LU_INDEX, lco, land_use, 0, 24, Land|use, 1, ranibow, $lco$, lco
    … …  
    285286lmaxth, lmaxth, upper_level_thermals, 0., 100., upper|level|thermals, 1, Greens, $top^{ĺev}_{th}$, toplev_th
    286287LLMAXTH, lmaxth, upper_level_thermals, 0., 100., upper|level|thermals, 1, Greens, $top^{ĺev}_{th}$, toplev_th
    287 lon, lon, longitude, -180., 180., longitude, degrees_East, seismic, $lon$, lon
    288 LON, lon, longitude, -180., 180., longitude, degrees_East, seismic, $lon$, lon
    289 west_east, lon, longitude, -180., 180., longitude, degrees_East, seismic, $lon$, lon
    290 XLONG, lon, longitude, -180., 180., longitude, degrees_East, seismic, $lon$, lon
    291 XLONG_M, lon, longitude, -180., 180., longitude, degrees_East, seismic, $lon$, lon
    292 longitude, lon, longitude, 0., 360., longitude, degrees_East, seismic, $lon$, lon
    293 lon, lon, longitude, 0., 360., longitude, degrees_East, seismic, $lon$, lon
    294 nav_lon, lon, longitude, 0., 360., longitude, degrees_East, seismic, $lon$, lon
     288lon, lon, longitude, -180., 180., longitude, degrees_east, seismic, $lon$, lon
     289LON, lon, longitude, -180., 180., longitude, degrees_east, seismic, $lon$, lon
     290west_east, lon, longitude, -180., 180., longitude, degrees_east, seismic, $lon$, lon
     291XLONG, lon, longitude, -180., 180., longitude, degrees_east, seismic, $lon$, lon
     292XLONG_M, lon, longitude, -180., 180., longitude, degrees_east, seismic, $lon$, lon
     293longitude, lon, longitude, 0., 360., longitude, degrees_east, seismic, $lon$, lon
     294lon, lon, longitude, 0., 360., longitude, degrees_east, seismic, $lon$, lon
     295nav_lon, lon, longitude, 0., 360., longitude, degrees_east, seismic, $lon$, lon
     296lon_bnds, lon_bnds, longitude_bnds, -180., 180., longitude_bounds, degrees_east, seismic, $lon_{bnds}$, lon_bnds
    295297mixthelaymean, mixthelaymean, mean_mixed_layer_potential_temperature, 0., 0.02, mean|mixed|layer|potential|temperature, K, Reds, $mixthelaymean$, mixthelaymean
    296298Mean mixed layer potential temperature, mixthelaymean, mean_mixed_layer_potential_temperature, 0., 0.02, mean|mixed|layer|potential|temperature, K, Reds, $mixthelaymean$, mixthelaymean
Note: See TracChangeset for help on using the changeset viewer.