Merge pull request #77 from wrf-model/release-v4.0.1

This merge adds fixes in the WPS v4.0.1 release onto the 'master' branch.
wrf-model · Oct 2, 2018 · 398505c · 398505c
2 parents 3ab4522 + 2266798
commit 398505c
Show file tree

Hide file tree

Showing 12 changed files with 126 additions and 74 deletions.
diff --git a/README b/README
@@ -1,4 +1,4 @@
-WRF Pre-Processing System Version 4.0 (8 June 2018)
+WRF Pre-Processing System Version 4.0.1 (2 October 2018)
 
 http://www2.mmm.ucar.edu/wrf/users/
 
@@ -10,7 +10,7 @@ programs that are part of WPS.  Both the ARW and NMM
 dynamical cores in WRF are supported by WPS.
 
 For questions and help to run the program, please see the 
-User's Guide at http://www2.mmm.ucar.edu/wrf/users/docs/user_guide_V3/contents.html
+User's Guide at http://www2.mmm.ucar.edu/wrf/users/docs/user_guide_v4/contents.html
 and send email to wrfhelp@ucar.edu.
 
 ===================================================

diff --git a/arch/preamble b/arch/preamble
@@ -33,14 +33,26 @@ PERL			=	perl
 RANLIB			=	echo
 
 #
-# Set path to compiled WRF code, which is assumed to be ../WRF for v4.x, or ../WRFV3 for v3.x
+# Look for compiled WRF code in one of several directories: WRF-4.0.1, WRF-4.0, WRF, and WRFV3.
+# The need for complicated logic to do this arises from the various names that the WRF code may take
+# on, depending on whether it was obtained through a GitHub archive file, as a clone of the GitHub
+# repository, or an older WRF v3.x tar file.
 # To override the path to the compiled WRF code, just set the WRF_DIR variable after the "endif" below
 #
-ifneq ($(wildcard $(DEV_TOP)/../WRF), ) # Check for WRF v4.x directory
-	WRF_DIR		=	../WRF
+ifneq ($(wildcard $(DEV_TOP)/../WRF-4.0.1), ) # Check for WRF v4.0.1 directory (probably GitHub archive)
+	WRF_DIR = ../WRF-4.0.1
 else
-	WRF_DIR		=	../WRFV3
+ifneq ($(wildcard $(DEV_TOP)/../WRF-4.0), )   # Check for WRF v4.0 directory (probably GitHub archive)
+	WRF_DIR = ../WRF-4.0
+else
+ifneq ($(wildcard $(DEV_TOP)/../WRF), )       # Check for clone of the WRF repository
+	WRF_DIR = ../WRF
+else                                          # Check for old WRF v3.x directory
+	WRF_DIR = ../WRFV3
+endif
 endif
+endif
+
 
 WRF_INCLUDE     =       -I$(WRF_DIR)/external/io_netcdf \
                         -I$(WRF_DIR)/external/io_grib_share \

diff --git a/geogrid/src/process_tile_module.F b/geogrid/src/process_tile_module.F
@@ -300,7 +300,7 @@ subroutine process_tile(which_domain, grid_type, dynopt,        &
       end if
 
       ! Initialize the output module now that we have the corner point lats/lons
-      call output_init(which_domain, 'OUTPUT FROM GEOGRID V4.0', '0000-00-00_00:00:00', grid_type, dynopt, &
+      call output_init(which_domain, 'OUTPUT FROM GEOGRID V4.0.1', '0000-00-00_00:00:00', grid_type, dynopt, &
                        corner_lats, corner_lons, &
                        start_dom_i,   end_dom_i,   start_dom_j,   end_dom_j, &
                        start_patch_i, end_patch_i, start_patch_j, end_patch_j, &

diff --git a/metgrid/src/Makefile b/metgrid/src/Makefile
@@ -66,7 +66,7 @@ queue_module.o: module_debug.o
 
 read_met_module.o: constants_module.o misc_definitions_module.o module_debug.o met_data_module.o
 
-rotate_winds_module.o: bitarray_module.o constants_module.o llxy_module.o misc_definitions_module.o module_debug.o
+rotate_winds_module.o: bitarray_module.o constants_module.o llxy_module.o misc_definitions_module.o module_debug.o module_map_utils.o
 
 storage_module.o: datatype_module.o minheap_module.o misc_definitions_module.o module_debug.o parallel_module.o module_stringutil.o
 

diff --git a/metgrid/src/input_module.F b/metgrid/src/input_module.F
@@ -427,7 +427,9 @@ subroutine read_global_attrs(title, start_date, grid_type, dyn_opt,
 #endif
 
          call ext_get_dom_ti_char          ('TITLE', title)
-         if (index(title,'GEOGRID V4.0') /= 0) then
+         if (index(title,'GEOGRID V4.0.1') /= 0) then
+            wps_version = 4.01
+         else if (index(title,'GEOGRID V4.0') /= 0) then
             wps_version = 4.0
          else if (index(title,'GEOGRID V3.9.1') /= 0) then
             wps_version = 3.91

diff --git a/metgrid/src/process_domain_module.F b/metgrid/src/process_domain_module.F
@@ -876,7 +876,7 @@ subroutine process_single_met_time(do_const_processing, &
       !   now we simply output every field from the storage module.
       !
 
-      title = 'OUTPUT FROM METGRID V4.0' 
+      title = 'OUTPUT FROM METGRID V4.0.1' 
 
       ! Initialize the output module for this domain and time
       call mprintf(.true.,LOGFILE,'Initializing output module.')

diff --git a/metgrid/src/rotate_winds_module.F b/metgrid/src/rotate_winds_module.F
@@ -3,6 +3,7 @@ module rotate_winds_module
    use bitarray_module
    use constants_module
    use llxy_module
+   use map_utils
    use misc_definitions_module
    use module_debug
 
@@ -98,6 +99,10 @@ subroutine metmap_xform(u, u_mask, v, v_mask, &
       real :: u_map, v_map, alpha, diff
       real, pointer, dimension(:,:) :: u_new, v_new, u_mult, v_mult
       logical :: do_last_col_u, do_last_row_u, do_last_col_v, do_last_row_v
+      real :: x_u, y_u, x_v, y_v
+      real :: xlat_u_p1 , xlon_u_p1, xlat_v_p1 , xlon_v_p1
+      real :: xlat_u_m1 , xlon_u_m1, xlat_v_m1 , xlon_v_m1
+      real :: diff_lon, diff_lat
 
       ! If the proj_info structure has not been initialized, we don't have
       !   information about the projection and standard longitude.
@@ -156,7 +161,6 @@ subroutine metmap_xform(u, u_mask, v, v_mask, &
             ! Rotate U field
             do j=us2,ue2
                do i=us1,ue1
-
                   diff = idir * (xlon_u(i,j) - proj_stack(current_nest_number)%stdlon)
                   if (diff > 180.) then
                      diff = diff - 360.
@@ -167,38 +171,55 @@ subroutine metmap_xform(u, u_mask, v, v_mask, &
                   ! Calculate the rotation angle, alpha, in radians
                   if (proj_stack(current_nest_number)%code == PROJ_LC) then
                      alpha = diff * proj_stack(current_nest_number)%cone * rad_per_deg * proj_stack(current_nest_number)%hemi 
-                  else if (proj_stack(current_nest_number)%code == PROJ_CASSINI) then
-                     if (j == ue2) then
-                        diff = xlon_u(i,j)-xlon_u(i,j-1)
-                        if (diff > 180.) then
-                           diff = diff - 360.
-                        else if (diff < -180.) then
-                           diff = diff + 360.
-                        end if
-                        alpha = atan2(   (-cos(xlat_u(i,j)*rad_per_deg) * diff*rad_per_deg),   &
-                                                        (xlat_u(i,j)-xlat_u(i,j-1))*rad_per_deg    &
-                                     )
-                     else if (j == us2) then
-                        diff = xlon_u(i,j+1)-xlon_u(i,j)
-                        if (diff > 180.) then
-                           diff = diff - 360.
-                        else if (diff < -180.) then
-                           diff = diff + 360.
+                  else if ( (proj_stack(SOURCE_PROJ)%code == PROJ_CASSINI) .and. ( idir == 1 ) ) then
+                        call latlon_to_ij ( proj_stack(SOURCE_PROJ) , &
+                                            xlat_u(i,j), xlon_u(i,j), x_u, y_u )
+                        call ij_to_latlon ( proj_stack(SOURCE_PROJ) , &
+                                            x_u, y_u + 0.1 , xlat_u_p1 , xlon_u_p1 )
+                        call ij_to_latlon ( proj_stack(SOURCE_PROJ) , &
+                                            x_u, y_u - 0.1 , xlat_u_m1 , xlon_u_m1 )
+                        diff_lon = xlon_u_p1-xlon_u_m1
+                        if (diff_lon > 180.) then
+                           diff_lon = diff_lon - 360.
+                        else if (diff_lon < -180.) then
+                           diff_lon = diff_lon + 360.
                         end if
-                        alpha = atan2(   (-cos(xlat_u(i,j)*rad_per_deg) * diff*rad_per_deg),   &
-                                                        (xlat_u(i,j+1)-xlat_u(i,j))*rad_per_deg    &
+                        diff_lat = xlat_u_p1-xlat_u_m1
+                        alpha =-atan2(   (-cos(xlat_u(i,j)*rad_per_deg) * diff_lon*rad_per_deg),   &
+                                                        diff_lat*rad_per_deg    &
                                      )
-                     else
-                        diff = xlon_u(i,j+1)-xlon_u(i,j-1)
-                        if (diff > 180.) then
-                           diff = diff - 360.
-                        else if (diff < -180.) then
-                           diff = diff + 360.
+                  else if ( (proj_stack(current_nest_number)%code == PROJ_CASSINI) .and. ( idir == -1 ) ) then
+                        if (j == ue2) then
+                           diff = xlon_u(i,j)-xlon_u(i,j-1)
+                           if (diff > 180.) then
+                              diff = diff - 360.
+                           else if (diff < -180.) then
+                              diff = diff + 360.
+                           end if
+                           alpha = atan2(   (-cos(xlat_u(i,j)*rad_per_deg) * diff*rad_per_deg),   &   
+                                                           (xlat_u(i,j)-xlat_u(i,j-1))*rad_per_deg    &   
+                                        )   
+                        else if (j == us2) then
+                           diff = xlon_u(i,j+1)-xlon_u(i,j)
+                           if (diff > 180.) then
+                              diff = diff - 360.
+                           else if (diff < -180.) then
+                              diff = diff + 360.
+                           end if
+                           alpha = atan2(   (-cos(xlat_u(i,j)*rad_per_deg) * diff*rad_per_deg),   &   
+                                                           (xlat_u(i,j+1)-xlat_u(i,j))*rad_per_deg    &   
+                                        )   
+                        else
+                           diff = xlon_u(i,j+1)-xlon_u(i,j-1)
+                           if (diff > 180.) then
+                              diff = diff - 360.
+                           else if (diff < -180.) then
+                              diff = diff + 360.
+                           end if
+                           alpha = atan2(   (-cos(xlat_u(i,j)*rad_per_deg) * diff*rad_per_deg),   &   
+                                                           (xlat_u(i,j+1)-xlat_u(i,j-1))*rad_per_deg    &   
+                                        )   
                         end if
-                        alpha = atan2(   (-cos(xlat_u(i,j)*rad_per_deg) * diff*rad_per_deg),   &
-                                                        (xlat_u(i,j+1)-xlat_u(i,j-1))*rad_per_deg    &
-                                     )
-                     end if
                   else
                      alpha = diff * rad_per_deg * proj_stack(current_nest_number)%hemi 
                   end if
@@ -257,38 +278,55 @@ subroutine metmap_xform(u, u_mask, v, v_mask, &
 
                   if (proj_stack(current_nest_number)%code == PROJ_LC) then
                      alpha = diff * proj_stack(current_nest_number)%cone * rad_per_deg * proj_stack(current_nest_number)%hemi 
-                  else if (proj_stack(current_nest_number)%code == PROJ_CASSINI) then
-                     if (j == ve2) then
-                        diff = xlon_v(i,j)-xlon_v(i,j-1)
-                        if (diff > 180.) then
-                           diff = diff - 360.
-                        else if (diff < -180.) then
-                           diff = diff + 360.
+                  else if ( (proj_stack(SOURCE_PROJ)%code == PROJ_CASSINI) .and. ( idir == 1 ) ) then
+                        call latlon_to_ij ( proj_stack(SOURCE_PROJ) , &
+                                            xlat_v(i,j), xlon_v(i,j), x_v, y_v )
+                        call ij_to_latlon ( proj_stack(SOURCE_PROJ) , &
+                                            x_v, y_v + 0.1 , xlat_v_p1 , xlon_v_p1 )
+                        call ij_to_latlon ( proj_stack(SOURCE_PROJ) , &
+                                            x_v, y_v - 0.1 , xlat_v_m1 , xlon_v_m1 )
+                        diff_lon = xlon_v_p1-xlon_v_m1
+                        if (diff_lon > 180.) then
+                           diff_lon = diff_lon - 360.
+                        else if (diff_lon < -180.) then
+                           diff_lon = diff_lon + 360.
                         end if
-                        alpha = atan2(   (-cos(xlat_v(i,j)*rad_per_deg) * diff*rad_per_deg),   &
-                                                        (xlat_v(i,j)-xlat_v(i,j-1))*rad_per_deg    &
+                        diff_lat = xlat_v_p1-xlat_v_m1
+                        alpha =-atan2(   (-cos(xlat_v(i,j)*rad_per_deg) * diff_lon*rad_per_deg),   &
+                                                        diff_lat*rad_per_deg    &
                                      )
-                     else if (j == vs2) then
-                        diff = xlon_v(i,j+1)-xlon_v(i,j)
-                        if (diff > 180.) then
-                           diff = diff - 360.
-                        else if (diff < -180.) then
-                           diff = diff + 360.
-                        end if
-                        alpha = atan2(   (-cos(xlat_v(i,j)*rad_per_deg) * diff*rad_per_deg),   &
-                                                        (xlat_v(i,j+1)-xlat_v(i,j))*rad_per_deg    &
-                                     )
-                     else
-                        diff = xlon_v(i,j+1)-xlon_v(i,j-1)
-                        if (diff > 180.) then
-                           diff = diff - 360.
-                        else if (diff < -180.) then
-                           diff = diff + 360.
+                  else if ( (proj_stack(current_nest_number)%code == PROJ_CASSINI) .and. ( idir == -1 ) ) then
+                        if (j == ve2) then
+                           diff = xlon_v(i,j)-xlon_v(i,j-1)
+                           if (diff > 180.) then
+                              diff = diff - 360.
+                           else if (diff < -180.) then
+                              diff = diff + 360.
+                           end if
+                           alpha = atan2(   (-cos(xlat_v(i,j)*rad_per_deg) * diff*rad_per_deg),   &
+                                                           (xlat_v(i,j)-xlat_v(i,j-1))*rad_per_deg    &
+                                        )
+                        else if (j == vs2) then
+                           diff = xlon_v(i,j+1)-xlon_v(i,j)
+                           if (diff > 180.) then
+                              diff = diff - 360.
+                           else if (diff < -180.) then
+                              diff = diff + 360.
+                           end if
+                           alpha = atan2(   (-cos(xlat_v(i,j)*rad_per_deg) * diff*rad_per_deg),   &
+                                                           (xlat_v(i,j+1)-xlat_v(i,j))*rad_per_deg    &
+                                        )
+                        else
+                           diff = xlon_v(i,j+1)-xlon_v(i,j-1)
+                           if (diff > 180.) then
+                              diff = diff - 360.
+                           else if (diff < -180.) then
+                              diff = diff + 360.
+                           end if
+                           alpha = atan2(   (-cos(xlat_v(i,j)*rad_per_deg) * diff*rad_per_deg),   &
+                                                           (xlat_v(i,j+1)-xlat_v(i,j-1))*rad_per_deg    &
+                                        )
                         end if
-                        alpha = atan2(   (-cos(xlat_v(i,j)*rad_per_deg) * diff*rad_per_deg),   &
-                                                        (xlat_v(i,j+1)-xlat_v(i,j-1))*rad_per_deg    &
-                                     )
-                     end if
                   else
                      alpha = diff * rad_per_deg * proj_stack(current_nest_number)%hemi 
                   end if

diff --git a/namelist.wps b/namelist.wps
@@ -36,7 +36,7 @@
  truelat1  =  30.0,
  truelat2  =  60.0,
  stand_lon = -98.0,
- geog_data_path = '/glade/p/work/wrfhelp/WPS_GEOG/'
+ geog_data_path = '/glade/work/wrfhelp/WPS_GEOG/'
 /
 
 &ungrib

diff --git a/namelist.wps.all_options b/namelist.wps.all_options
@@ -59,7 +59,7 @@
  truelat1  =  30.0,
  truelat2  =  60.0,
  stand_lon = -98.0,
- geog_data_path = '/glade/p/work/wrfhelp/WPS_GEOG/'
+ geog_data_path = '/glade/work/wrfhelp/WPS_GEOG/'
  opt_geogrid_tbl_path = 'geogrid/'
 /
  geog_data_res     = 'modis_lakes+10m','modis_lakes+2m',

diff --git a/namelist.wps.fire b/namelist.wps.fire
@@ -42,7 +42,7 @@
  truelat1  =  39.338,
  truelat2  =  39.338,
  stand_lon = -106.807,
- geog_data_path = '/glade/p/work/wrfhelp/WPS_GEOG/'
+ geog_data_path = '/glade/work/wrfhelp/WPS_GEOG/'
 /
 
 &ungrib

diff --git a/namelist.wps.global b/namelist.wps.global
@@ -32,7 +32,7 @@
  stand_lon = 0.
  pole_lat = 90.0
  pole_lon = 0.0
- geog_data_path = '/glade/p/work/wrfhelp/WPS_GEOG/'
+ geog_data_path = '/glade/work/wrfhelp/WPS_GEOG/'
 /
  ref_lat = 45.0
  ref_lon = -98

diff --git a/namelist.wps.nmm b/namelist.wps.nmm
@@ -20,7 +20,7 @@
  map_proj = 'rotated_ll',
  ref_lat   =  42.00,
  ref_lon   = -71.00,
- geog_data_path = '/glade/p/work/wrfhelp/WPS_GEOG/'
+ geog_data_path = '/glade/work/wrfhelp/WPS_GEOG/'
 /
 
 &ungrib