Changeset 3459 for palm/trunk/SOURCE
- Timestamp:
- Oct 30, 2018 3:04:11 PM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 2 edited
-
netcdf_data_input_mod.f90 (modified) (7 diffs)
-
netcdf_interface_mod.f90 (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/netcdf_data_input_mod.f90
r3458 r3459 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Define coordinate reference system (crs) and read from input dataset 23 ! Revise default values for reference coordinates 23 24 ! 24 25 ! Former revisions: … … 348 349 REAL(wp) :: fill_v !< fill value for v 349 350 REAL(wp) :: fill_w !< fill value for w 350 REAL(wp) :: latitude !< latitude of the southern model boundary351 REAL(wp) :: longitude !< longitude of the western model boundary352 REAL(wp) :: origin_x = 0.0_wp !< x position of the western model boundary (default: Hannover)353 REAL(wp) :: origin_y = 0.0_wp !< y position of the northern model boundary (default: Hannover)354 REAL(wp) :: origin_z = 0.0_wp !< reference height of input data (default: Hannover)351 REAL(wp) :: latitude = 0.0_wp !< latitude of the lower left corner 352 REAL(wp) :: longitude = -3.0_wp !< longitude of the lower left corner 353 REAL(wp) :: origin_x = 500000.0_wp !< UTM easting of the lower left corner 354 REAL(wp) :: origin_y = 0.0_wp !< UTM northing of the lower left corner 355 REAL(wp) :: origin_z = 0.0_wp !< reference height of input data 355 356 REAL(wp) :: rotation_angle = 0.0_wp !< rotation angle of input data 356 357 … … 585 586 END TYPE global_atts_type 586 587 ! 588 !-- Define type for coordinate reference system (crs) 589 TYPE crs_type 590 CHARACTER(LEN=200) :: epsg_code = 'EPSG:25831' !< EPSG code 591 CHARACTER(LEN=200) :: grid_mapping_name = 'transverse_mercator' !< name of grid mapping 592 CHARACTER(LEN=200) :: long_name = 'coordinate reference system' !< name of variable crs 593 CHARACTER(LEN=200) :: units = 'm' !< unit of crs 594 595 REAL(wp) :: false_easting = 500000.0_wp !< false easting 596 REAL(wp) :: false_northing = 0.0_wp !< false northing 597 REAL(wp) :: inverse_flattening = 298.257223563_wp !< 1/f (default for WGS84) 598 REAL(wp) :: latitude_of_projection_origin = 0.0_wp !< latitude of projection origin 599 REAL(wp) :: longitude_of_central_meridian = 3.0_wp !< longitude of central meridian of UTM zone (default: zone 31) 600 REAL(wp) :: longitude_of_prime_meridian = 0.0_wp !< longitude of prime meridian 601 REAL(wp) :: scale_factor_at_central_meridian = 0.9996_wp !< scale factor of UTM coordinates 602 REAL(wp) :: semi_major_axis = 6378137.0_wp !< length of semi major axis (default for WGS84) 603 END TYPE crs_type 604 605 ! 587 606 !-- Define variables 588 TYPE(dims_xy) :: dim_static !< data structure for x, y-dimension in static input file 607 TYPE(crs_type) :: coord_ref_sys !< coordinate reference system 608 609 TYPE(dims_xy) :: dim_static !< data structure for x, y-dimension in static input file 589 610 590 611 TYPE(nest_offl_type) :: nest_offl !< data structure for data input at lateral and top boundaries (provided by Inifor) … … 740 761 building_id_f, building_pars_f, building_type_f, & 741 762 chem_emis, chem_emis_att, chem_emis_att_type, chem_emis_val_type, & 763 coord_ref_sys, & 742 764 init_3d, init_model, input_file_static, input_pids_static, & 743 765 input_pids_dynamic, leaf_area_density_f, nest_offl, & … … 796 818 ! Description: 797 819 ! ------------ 798 !> Reads global attributes required for initialization of the model. 820 !> Reads global attributes and coordinate reference system required for 821 !> initialization of the model. 799 822 !------------------------------------------------------------------------------! 800 823 SUBROUTINE netcdf_data_input_init … … 802 825 IMPLICIT NONE 803 826 804 INTEGER(iwp) :: id_mod !< NetCDF id of input file 827 INTEGER(iwp) :: id_mod !< NetCDF id of input file 828 INTEGER(iwp) :: var_id_crs !< NetCDF id of variable crs 805 829 806 830 IF ( .NOT. input_pids_static ) RETURN … … 833 857 CALL get_attribute( id_mod, input_file_atts%rotation_angle_char, & 834 858 input_file_atts%rotation_angle, .TRUE. ) 835 859 ! 860 !-- Read coordinate reference system if available 861 nc_stat = NF90_INQ_VARID( id_mod, 'crs', var_id_crs ) 862 IF ( nc_stat == NF90_NOERR ) THEN 863 CALL get_attribute( id_mod, 'epsg_code', & 864 coord_ref_sys%epsg_code, & 865 .FALSE., 'crs' ) 866 CALL get_attribute( id_mod, 'false_easting', & 867 coord_ref_sys%false_easting, & 868 .FALSE., 'crs' ) 869 CALL get_attribute( id_mod, 'false_northing', & 870 coord_ref_sys%false_northing, & 871 .FALSE., 'crs' ) 872 CALL get_attribute( id_mod, 'grid_mapping_name', & 873 coord_ref_sys%grid_mapping_name, & 874 .FALSE., 'crs' ) 875 CALL get_attribute( id_mod, 'inverse_flattening', & 876 coord_ref_sys%inverse_flattening, & 877 .FALSE., 'crs' ) 878 CALL get_attribute( id_mod, 'latitude_of_projection_origin', & 879 coord_ref_sys%latitude_of_projection_origin, & 880 .FALSE., 'crs' ) 881 CALL get_attribute( id_mod, 'long_name', & 882 coord_ref_sys%long_name, & 883 .FALSE., 'crs' ) 884 CALL get_attribute( id_mod, 'longitude_of_central_meridian', & 885 coord_ref_sys%longitude_of_central_meridian, & 886 .FALSE., 'crs' ) 887 CALL get_attribute( id_mod, 'longitude_of_prime_meridian', & 888 coord_ref_sys%longitude_of_prime_meridian, & 889 .FALSE., 'crs' ) 890 CALL get_attribute( id_mod, 'scale_factor_at_central_meridian', & 891 coord_ref_sys%scale_factor_at_central_meridian, & 892 .FALSE., 'crs' ) 893 CALL get_attribute( id_mod, 'semi_major_axis', & 894 coord_ref_sys%semi_major_axis, & 895 .FALSE., 'crs' ) 896 CALL get_attribute( id_mod, 'units', & 897 coord_ref_sys%units, & 898 .FALSE., 'crs' ) 899 ENDIF 836 900 ! 837 901 !-- Finally, close input file -
palm/trunk/SOURCE/netcdf_interface_mod.f90
r3448 r3459 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! - Add variable crs to output files 23 ! - Add long_name to UTM coordinates 24 ! - Add latitude/longitude coordinates. For 3d and xy output, lon and lat are 25 ! only written if parallel output is used. 23 26 ! 24 27 ! Former revisions: … … 311 314 !> @todo timeseries and profile output still needs to be rewritten to allow 312 315 !> modularization 316 !> @todo output 2d UTM coordinates without global arrays 317 !> @todo output longitude/latitude also with non-parallel output (3d and xy) 313 318 !------------------------------------------------------------------------------! 314 319 MODULE netcdf_interface … … 324 329 USE mas_global_attributes, & 325 330 ONLY: dim_size_agtnum 331 332 USE netcdf_data_input_mod, & 333 ONLY: coord_ref_sys 326 334 327 335 PRIVATE … … 472 480 id_var_eutm_xy, id_var_nutm_xy, & 473 481 id_var_eutm_xz, id_var_nutm_xz, & 474 id_var_eutm_yz, id_var_nutm_yz, & 475 id_var_eutm_ma, id_var_nutm_ma 482 id_var_eutm_yz, id_var_nutm_yz 483 484 INTEGER(iwp), DIMENSION(0:1,0:2) :: id_var_lat_3d, id_var_lon_3d, & 485 id_var_lat_xy, id_var_lon_xy, & 486 id_var_lat_xz, id_var_lon_xz, & 487 id_var_lat_yz, id_var_lon_yz 476 488 477 489 INTEGER :: netcdf_data_format = 2 !< NetCDF3 64bit offset format … … 520 532 INTEGER(iwp), DIMENSION(1:max_masks,0:1,0:2) :: id_var_eutm_mask, & 521 533 id_var_nutm_mask 534 535 INTEGER(iwp), DIMENSION(1:max_masks,0:1,0:2) :: id_var_lat_mask, & 536 id_var_lon_mask 522 537 523 538 INTEGER(iwp), DIMENSION(1:max_masks,0:1,100) :: id_var_domask … … 624 639 625 640 USE netcdf_data_input_mod, & 626 ONLY: init_model641 ONLY: init_model 627 642 628 643 USE ocean_mod, & … … 720 735 721 736 REAL(wp) :: cos_ra !< cosine of rotation_angle 737 REAL(wp) :: eutm !< easting (UTM) 738 REAL(wp) :: nutm !< northing (UTM) 722 739 REAL(wp) :: shift_x !< shift of x coordinate 723 740 REAL(wp) :: shift_y !< shift of y coordinate … … 725 742 726 743 REAL(wp), DIMENSION(1) :: last_time_coordinate !< last time value in file 744 REAL(wp), DIMENSION(8) :: crs_list !< list of coord_ref_sys values 727 745 728 746 REAL(wp), DIMENSION(:), ALLOCATABLE :: netcdf_data !< 729 747 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: netcdf_data_2d !< 748 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lat !< latitude 749 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lon !< longitude 730 750 731 751 … … 822 842 823 843 ! 844 !-- Convert coord_ref_sys into vector (used for lat/lon calculation) 845 crs_list = (/ coord_ref_sys%semi_major_axis, & 846 coord_ref_sys%inverse_flattening, & 847 coord_ref_sys%longitude_of_prime_meridian, & 848 coord_ref_sys%longitude_of_central_meridian, & 849 coord_ref_sys%scale_factor_at_central_meridian, & 850 coord_ref_sys%latitude_of_projection_origin, & 851 coord_ref_sys%false_easting, & 852 coord_ref_sys%false_northing /) 853 854 ! 824 855 !-- Determine the mode to be processed 825 856 IF ( extend ) THEN … … 962 993 (/ id_dim_x_mask(mid,av) /), & 963 994 'E_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,0), & 964 'm', ' ', 000, 000, 000 )995 'm', 'projection_x_coordinate', 000, 000, 000 ) 965 996 CALL netcdf_create_var( id_set_mask(mid,av), & 966 997 (/ id_dim_y_mask(mid,av) /), & 967 998 'N_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,0), & 968 'm', ' ', 000, 000, 000 )999 'm', 'projection_y_coordinate', 000, 000, 000 ) 969 1000 CALL netcdf_create_var( id_set_mask(mid,av), & 970 1001 (/ id_dim_xu_mask(mid,av) /), & 971 1002 'Eu_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,1), & 972 'm', ' ', 000, 000, 000 )1003 'm', 'projection_x_coordinate', 000, 000, 000 ) 973 1004 CALL netcdf_create_var( id_set_mask(mid,av), & 974 1005 (/ id_dim_y_mask(mid,av) /), & 975 1006 'Nu_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,1), & 976 'm', ' ', 000, 000, 000 )1007 'm', 'projection_y_coordinate', 000, 000, 000 ) 977 1008 CALL netcdf_create_var( id_set_mask(mid,av), & 978 1009 (/ id_dim_x_mask(mid,av) /), & 979 1010 'Ev_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,2), & 980 'm', ' ', 000, 000, 000 )1011 'm', 'projection_x_coordinate', 000, 000, 000 ) 981 1012 CALL netcdf_create_var( id_set_mask(mid,av), & 982 1013 (/ id_dim_yv_mask(mid,av) /), & 983 1014 'Nv_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,2), & 984 'm', ' ', 000, 000, 000 )1015 'm', 'projection_y_coordinate', 000, 000, 000 ) 985 1016 ELSE 986 1017 CALL netcdf_create_var( id_set_mask(mid,av), & 987 (/ id_dim_x_mask(mid,av), id_dim_y_mask(mid,av) /), 988 'E_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,0), &989 'm', ' ', 000, 000, 000 )1018 (/ id_dim_x_mask(mid,av), id_dim_y_mask(mid,av) /), & 1019 'E_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,0), & 1020 'm', 'projection_x_coordinate', 000, 000, 000 ) 990 1021 CALL netcdf_create_var( id_set_mask(mid,av), & 991 (/ id_dim_x_mask(mid,av), id_dim_y_mask(mid,av) /), 992 'N_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,0), &993 'm', ' ', 000, 000, 000 )1022 (/ id_dim_x_mask(mid,av), id_dim_y_mask(mid,av) /), & 1023 'N_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,0), & 1024 'm', 'projection_y_coordinate', 000, 000, 000 ) 994 1025 CALL netcdf_create_var( id_set_mask(mid,av), & 995 (/ id_dim_xu_mask(mid,av), id_dim_y_mask(mid,av) /), 996 'Eu_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,1), &997 'm', ' ', 000, 000, 000 )1026 (/ id_dim_xu_mask(mid,av), id_dim_y_mask(mid,av) /),& 1027 'Eu_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,1), & 1028 'm', 'projection_x_coordinate', 000, 000, 000 ) 998 1029 CALL netcdf_create_var( id_set_mask(mid,av), & 999 (/ id_dim_xu_mask(mid,av), id_dim_y_mask(mid,av) /), 1000 'Nu_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,1), &1001 'm', ' ', 000, 000, 000 )1030 (/ id_dim_xu_mask(mid,av), id_dim_y_mask(mid,av) /),& 1031 'Nu_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,1), & 1032 'm', 'projection_y_coordinate', 000, 000, 000 ) 1002 1033 CALL netcdf_create_var( id_set_mask(mid,av), & 1003 (/ id_dim_x_mask(mid,av), id_dim_yv_mask(mid,av) /), 1004 'Ev_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,2), &1005 'm', ' ', 000, 000, 000 )1034 (/ id_dim_x_mask(mid,av), id_dim_yv_mask(mid,av) /),& 1035 'Ev_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,2), & 1036 'm', 'projection_x_coordinate', 000, 000, 000 ) 1006 1037 CALL netcdf_create_var( id_set_mask(mid,av), & 1007 (/ id_dim_x_mask(mid,av), id_dim_yv_mask(mid,av) /), & 1008 'Nv_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,2), & 1009 'm', '', 000, 000, 000 ) 1010 ENDIF 1038 (/ id_dim_x_mask(mid,av), id_dim_yv_mask(mid,av) /),& 1039 'Nv_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,2), & 1040 'm', 'projection_y_coordinate', 000, 000, 000 ) 1041 ENDIF 1042 ! 1043 !-- Define geographic coordinates 1044 CALL netcdf_create_var( id_set_mask(mid,av), & 1045 (/ id_dim_x_mask(mid,av), id_dim_y_mask(mid,av) /), & 1046 'lon', NF90_DOUBLE, id_var_lon_mask(mid,av,0), & 1047 'degrees_east', 'longitude', 000, 000, 000 ) 1048 CALL netcdf_create_var( id_set_mask(mid,av), & 1049 (/ id_dim_x_mask(mid,av), id_dim_y_mask(mid,av) /), & 1050 'lat', NF90_DOUBLE, id_var_lat_mask(mid,av,0), & 1051 'degrees_north', 'latitude', 000, 000, 000 ) 1052 CALL netcdf_create_var( id_set_mask(mid,av), & 1053 (/ id_dim_xu_mask(mid,av), id_dim_y_mask(mid,av) /),& 1054 'lonu', NF90_DOUBLE, id_var_lon_mask(mid,av,1), & 1055 'degrees_east', 'longitude', 000, 000, 000 ) 1056 CALL netcdf_create_var( id_set_mask(mid,av), & 1057 (/ id_dim_xu_mask(mid,av), id_dim_y_mask(mid,av) /),& 1058 'latu', NF90_DOUBLE, id_var_lat_mask(mid,av,1), & 1059 'degrees_north', 'latitude', 000, 000, 000 ) 1060 CALL netcdf_create_var( id_set_mask(mid,av), & 1061 (/ id_dim_x_mask(mid,av), id_dim_yv_mask(mid,av) /),& 1062 'lonv', NF90_DOUBLE, id_var_lon_mask(mid,av,2), & 1063 'degrees_east', 'longitude', 000, 000, 000 ) 1064 CALL netcdf_create_var( id_set_mask(mid,av), & 1065 (/ id_dim_x_mask(mid,av), id_dim_yv_mask(mid,av) /),& 1066 'latv', NF90_DOUBLE, id_var_lat_mask(mid,av,2), & 1067 'degrees_north', 'latitude', 000, 000, 000 ) 1068 ! 1069 !-- Define coordinate-reference system 1070 CALL netcdf_create_crs( id_set_mask(mid,av), 000 ) 1011 1071 ! 1012 1072 !-- In case of non-flat topography define 2d-arrays containing the height … … 1364 1424 DEALLOCATE( netcdf_data_2d ) 1365 1425 ENDIF 1366 1426 ! 1427 !-- Write lon and lat data. 1428 ALLOCATE( lat(mask_size(mid,1),mask_size(mid,2)) ) 1429 ALLOCATE( lon(mask_size(mid,1),mask_size(mid,2)) ) 1430 cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp ) 1431 sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp ) 1432 1433 DO k = 0, 2 1434 ! 1435 !-- Scalar grid points 1436 IF ( k == 0 ) THEN 1437 shift_x = 0.5 ; shift_y = 0.5 1438 ! 1439 !-- u grid points 1440 ELSEIF ( k == 1 ) THEN 1441 shift_x = 0.0 ; shift_y = 0.5 1442 ! 1443 !-- v grid points 1444 ELSEIF ( k == 2 ) THEN 1445 shift_x = 0.5 ; shift_y = 0.0 1446 ENDIF 1447 1448 DO j = 1, mask_size(mid,2) 1449 DO i = 1, mask_size(mid,1) 1450 eutm = init_model%origin_x & 1451 + cos_ra * ( mask_i_global(mid,i) + shift_x ) * dx & 1452 + sin_ra * ( mask_j_global(mid,j) + shift_y ) * dy 1453 nutm = init_model%origin_y & 1454 - sin_ra * ( mask_i_global(mid,i) + shift_x ) * dx & 1455 + cos_ra * ( mask_j_global(mid,j) + shift_y ) * dy 1456 1457 CALL convert_utm_to_geographic( crs_list, & 1458 eutm, nutm, & 1459 lon(i,j), lat(i,j) ) 1460 ENDDO 1461 ENDDO 1462 1463 nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), & 1464 id_var_lon_mask(mid,av,k), & 1465 lon, start = (/ 1, 1 /), & 1466 count = (/ mask_size(mid,1), & 1467 mask_size(mid,2) /) ) 1468 CALL netcdf_handle_error( 'netcdf_define_header', 556 ) 1469 1470 nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), & 1471 id_var_lat_mask(mid,av,k), & 1472 lat, start = (/ 1, 1 /), & 1473 count = (/ mask_size(mid,1), & 1474 mask_size(mid,2) /) ) 1475 CALL netcdf_handle_error( 'netcdf_define_header', 556 ) 1476 ENDDO 1477 1478 DEALLOCATE( lat ) 1479 DEALLOCATE( lon ) 1367 1480 ! 1368 1481 !-- Write zu and zw data (vertical axes) … … 1703 1816 (/ id_dim_x_3d(av) /), & 1704 1817 'E_UTM', NF90_DOUBLE, id_var_eutm_3d(av,0), & 1705 'm', ' ', 000, 000, 000 )1818 'm', 'projection_x_coordinate', 000, 000, 000 ) 1706 1819 CALL netcdf_create_var( id_set_3d(av), & 1707 1820 (/ id_dim_y_3d(av) /), & 1708 1821 'N_UTM', NF90_DOUBLE, id_var_nutm_3d(av,0), & 1709 'm', ' ', 000, 000, 000 )1822 'm', 'projection_y_coordinate', 000, 000, 000 ) 1710 1823 CALL netcdf_create_var( id_set_3d(av), & 1711 1824 (/ id_dim_xu_3d(av) /), & 1712 1825 'Eu_UTM', NF90_DOUBLE, id_var_eutm_3d(av,1), & 1713 'm', ' ', 000, 000, 000 )1826 'm', 'projection_x_coordinate', 000, 000, 000 ) 1714 1827 CALL netcdf_create_var( id_set_3d(av), & 1715 1828 (/ id_dim_y_3d(av) /), & 1716 1829 'Nu_UTM', NF90_DOUBLE, id_var_nutm_3d(av,1), & 1717 'm', ' ', 000, 000, 000 )1830 'm', 'projection_y_coordinate', 000, 000, 000 ) 1718 1831 CALL netcdf_create_var( id_set_3d(av), & 1719 1832 (/ id_dim_x_3d(av) /), & 1720 1833 'Ev_UTM', NF90_DOUBLE, id_var_eutm_3d(av,2), & 1721 'm', ' ', 000, 000, 000 )1834 'm', 'projection_x_coordinate', 000, 000, 000 ) 1722 1835 CALL netcdf_create_var( id_set_3d(av), & 1723 1836 (/ id_dim_yv_3d(av) /), & 1724 1837 'Nv_UTM', NF90_DOUBLE, id_var_nutm_3d(av,2), & 1725 'm', ' ', 000, 000, 000 )1838 'm', 'projection_y_coordinate', 000, 000, 000 ) 1726 1839 ELSE 1727 1840 CALL netcdf_create_var( id_set_3d(av), & 1728 1841 (/ id_dim_x_3d(av), id_dim_y_3d(av) /), & 1729 1842 'E_UTM', NF90_DOUBLE, id_var_eutm_3d(av,0), & 1730 'm', ' ', 000, 000, 000 )1843 'm', 'projection_x_coordinate', 000, 000, 000 ) 1731 1844 CALL netcdf_create_var( id_set_3d(av), & 1732 1845 (/ id_dim_x_3d(av), id_dim_y_3d(av) /), & 1733 1846 'N_UTM', NF90_DOUBLE, id_var_nutm_3d(av,0), & 1734 'm', ' ', 000, 000, 000 )1847 'm', 'projection_y_coordinate', 000, 000, 000 ) 1735 1848 CALL netcdf_create_var( id_set_3d(av), & 1736 1849 (/ id_dim_xu_3d(av), id_dim_y_3d(av) /), & 1737 1850 'Eu_UTM', NF90_DOUBLE, id_var_eutm_3d(av,1), & 1738 'm', ' ', 000, 000, 000 )1851 'm', 'projection_x_coordinate', 000, 000, 000 ) 1739 1852 CALL netcdf_create_var( id_set_3d(av), & 1740 1853 (/ id_dim_xu_3d(av), id_dim_y_3d(av) /), & 1741 1854 'Nu_UTM', NF90_DOUBLE, id_var_nutm_3d(av,1), & 1742 'm', ' ', 000, 000, 000 )1855 'm', 'projection_y_coordinate', 000, 000, 000 ) 1743 1856 CALL netcdf_create_var( id_set_3d(av), & 1744 1857 (/ id_dim_x_3d(av), id_dim_yv_3d(av) /), & 1745 1858 'Ev_UTM', NF90_DOUBLE, id_var_eutm_3d(av,2), & 1746 'm', ' ', 000, 000, 000 )1859 'm', 'projection_x_coordinate', 000, 000, 000 ) 1747 1860 CALL netcdf_create_var( id_set_3d(av), & 1748 1861 (/ id_dim_x_3d(av), id_dim_yv_3d(av) /), & 1749 1862 'Nv_UTM', NF90_DOUBLE, id_var_nutm_3d(av,2), & 1750 'm', '', 000, 000, 000 ) 1751 ENDIF 1863 'm', 'projection_y_coordinate', 000, 000, 000 ) 1864 ENDIF 1865 ! 1866 !-- Define geographic coordinates 1867 CALL netcdf_create_var( id_set_3d(av), & 1868 (/ id_dim_x_3d(av), id_dim_y_3d(av) /), & 1869 'lon', NF90_DOUBLE, id_var_lon_3d(av,0), & 1870 'degrees_east', 'longitude', 000, 000, 000 ) 1871 CALL netcdf_create_var( id_set_3d(av), & 1872 (/ id_dim_x_3d(av), id_dim_y_3d(av) /), & 1873 'lat', NF90_DOUBLE, id_var_lat_3d(av,0), & 1874 'degrees_north', 'latitude', 000, 000, 000 ) 1875 CALL netcdf_create_var( id_set_3d(av), & 1876 (/ id_dim_xu_3d(av), id_dim_y_3d(av) /), & 1877 'lonu', NF90_DOUBLE, id_var_lon_3d(av,1), & 1878 'degrees_east', 'longitude', 000, 000, 000 ) 1879 CALL netcdf_create_var( id_set_3d(av), & 1880 (/ id_dim_xu_3d(av), id_dim_y_3d(av) /), & 1881 'latu', NF90_DOUBLE, id_var_lat_3d(av,1), & 1882 'degrees_north', 'latitude', 000, 000, 000 ) 1883 CALL netcdf_create_var( id_set_3d(av), & 1884 (/ id_dim_x_3d(av), id_dim_yv_3d(av) /), & 1885 'lonv', NF90_DOUBLE, id_var_lon_3d(av,2), & 1886 'degrees_east', 'longitude', 000, 000, 000 ) 1887 CALL netcdf_create_var( id_set_3d(av), & 1888 (/ id_dim_x_3d(av), id_dim_yv_3d(av) /), & 1889 'latv', NF90_DOUBLE, id_var_lat_3d(av,2), & 1890 'degrees_north', 'latitude', 000, 000, 000 ) 1891 ! 1892 !-- Define coordinate-reference system 1893 CALL netcdf_create_crs( id_set_3d(av), 000 ) 1752 1894 ! 1753 1895 !-- In case of non-flat topography define 2d-arrays containing the height … … 2168 2310 CALL netcdf_handle_error( 'netcdf_define_header', 86 ) 2169 2311 ENDIF 2312 2313 ENDIF 2314 ! 2315 !-- Write lon and lat data. Only for parallel output. 2316 IF ( netcdf_data_format > 4 ) THEN 2317 2318 ALLOCATE( lat(nxl:nxr,nys:nyn) ) 2319 ALLOCATE( lon(nxl:nxr,nys:nyn) ) 2320 cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp ) 2321 sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp ) 2322 2323 DO k = 0, 2 2324 ! 2325 !-- Scalar grid points 2326 IF ( k == 0 ) THEN 2327 shift_x = 0.5 ; shift_y = 0.5 2328 ! 2329 !-- u grid points 2330 ELSEIF ( k == 1 ) THEN 2331 shift_x = 0.0 ; shift_y = 0.5 2332 ! 2333 !-- v grid points 2334 ELSEIF ( k == 2 ) THEN 2335 shift_x = 0.5 ; shift_y = 0.0 2336 ENDIF 2337 2338 DO j = nys, nyn 2339 DO i = nxl, nxr 2340 eutm = init_model%origin_x & 2341 + cos_ra * ( i + shift_x ) * dx & 2342 + sin_ra * ( j + shift_y ) * dy 2343 nutm = init_model%origin_y & 2344 - sin_ra * ( i + shift_x ) * dx & 2345 + cos_ra * ( j + shift_y ) * dy 2346 2347 CALL convert_utm_to_geographic( crs_list, & 2348 eutm, nutm, & 2349 lon(i,j), lat(i,j) ) 2350 ENDDO 2351 ENDDO 2352 2353 nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_lon_3d(av,k), & 2354 lon, start = (/ nxl+1, nys+1 /), & 2355 count = (/ nxr-nxl+1, nyn-nys+1 /) ) 2356 CALL netcdf_handle_error( 'netcdf_define_header', 556 ) 2357 2358 nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_lat_3d(av,k), & 2359 lat, start = (/ nxl+1, nys+1 /), & 2360 count = (/ nxr-nxl+1, nyn-nys+1 /) ) 2361 CALL netcdf_handle_error( 'netcdf_define_header', 556 ) 2362 ENDDO 2363 2364 DEALLOCATE( lat ) 2365 DEALLOCATE( lon ) 2170 2366 2171 2367 ENDIF … … 2673 2869 (/ id_dim_x_xy(av) /), & 2674 2870 'E_UTM', NF90_DOUBLE, id_var_eutm_xy(av,0), & 2675 'm', ' ', 000, 000, 000 )2871 'm', 'projection_x_coordinate', 000, 000, 000 ) 2676 2872 CALL netcdf_create_var( id_set_xy(av), & 2677 2873 (/ id_dim_y_xy(av) /), & 2678 2874 'N_UTM', NF90_DOUBLE, id_var_nutm_xy(av,0), & 2679 'm', ' ', 000, 000, 000 )2875 'm', 'projection_y_coordinate', 000, 000, 000 ) 2680 2876 CALL netcdf_create_var( id_set_xy(av), & 2681 2877 (/ id_dim_xu_xy(av) /), & 2682 2878 'Eu_UTM', NF90_DOUBLE, id_var_eutm_xy(av,1), & 2683 'm', ' ', 000, 000, 000 )2879 'm', 'projection_x_coordinate', 000, 000, 000 ) 2684 2880 CALL netcdf_create_var( id_set_xy(av), & 2685 2881 (/ id_dim_y_xy(av) /), & 2686 2882 'Nu_UTM', NF90_DOUBLE, id_var_nutm_xy(av,1), & 2687 'm', ' ', 000, 000, 000 )2883 'm', 'projection_y_coordinate', 000, 000, 000 ) 2688 2884 CALL netcdf_create_var( id_set_xy(av), & 2689 2885 (/ id_dim_x_xy(av) /), & 2690 2886 'Ev_UTM', NF90_DOUBLE, id_var_eutm_xy(av,2), & 2691 'm', ' ', 000, 000, 000 )2887 'm', 'projection_x_coordinate', 000, 000, 000 ) 2692 2888 CALL netcdf_create_var( id_set_xy(av), & 2693 2889 (/ id_dim_yv_xy(av) /), & 2694 2890 'Nv_UTM', NF90_DOUBLE, id_var_nutm_xy(av,2), & 2695 'm', ' ', 000, 000, 000 )2891 'm', 'projection_y_coordinate', 000, 000, 000 ) 2696 2892 ELSE 2697 2893 CALL netcdf_create_var( id_set_xy(av), & 2698 2894 (/ id_dim_x_xy(av), id_dim_y_xy(av) /), & 2699 2895 'E_UTM', NF90_DOUBLE, id_var_eutm_xy(av,0), & 2700 'm', ' ', 000, 000, 000 )2896 'm', 'projection_x_coordinate', 000, 000, 000 ) 2701 2897 CALL netcdf_create_var( id_set_xy(av), & 2702 2898 (/ id_dim_x_xy(av), id_dim_y_xy(av) /), & 2703 2899 'N_UTM', NF90_DOUBLE, id_var_nutm_xy(av,0), & 2704 'm', ' ', 000, 000, 000 )2900 'm', 'projection_y_coordinate', 000, 000, 000 ) 2705 2901 CALL netcdf_create_var( id_set_xy(av), & 2706 2902 (/ id_dim_xu_xy(av), id_dim_y_xy(av) /), & 2707 2903 'Eu_UTM', NF90_DOUBLE, id_var_eutm_xy(av,1), & 2708 'm', ' ', 000, 000, 000 )2904 'm', 'projection_x_coordinate', 000, 000, 000 ) 2709 2905 CALL netcdf_create_var( id_set_xy(av), & 2710 2906 (/ id_dim_xu_xy(av), id_dim_y_xy(av) /), & 2711 2907 'Nu_UTM', NF90_DOUBLE, id_var_nutm_xy(av,1), & 2712 'm', ' ', 000, 000, 000 )2908 'm', 'projection_y_coordinate', 000, 000, 000 ) 2713 2909 CALL netcdf_create_var( id_set_xy(av), & 2714 2910 (/ id_dim_x_xy(av), id_dim_yv_xy(av) /), & 2715 2911 'Ev_UTM', NF90_DOUBLE, id_var_eutm_xy(av,2), & 2716 'm', ' ', 000, 000, 000 )2912 'm', 'projection_x_coordinate', 000, 000, 000 ) 2717 2913 CALL netcdf_create_var( id_set_xy(av), & 2718 2914 (/ id_dim_x_xy(av), id_dim_yv_xy(av) /), & 2719 2915 'Nv_UTM', NF90_DOUBLE, id_var_nutm_xy(av,2), & 2720 'm', '', 000, 000, 000 ) 2721 ENDIF 2916 'm', 'projection_y_coordinate', 000, 000, 000 ) 2917 ENDIF 2918 ! 2919 !-- Define geographic coordinates 2920 CALL netcdf_create_var( id_set_xy(av), & 2921 (/ id_dim_x_xy(av), id_dim_y_xy(av) /), & 2922 'lon', NF90_DOUBLE, id_var_lon_xy(av,0), & 2923 'degrees_east', 'longitude', 000, 000, 000 ) 2924 CALL netcdf_create_var( id_set_xy(av), & 2925 (/ id_dim_x_xy(av), id_dim_y_xy(av) /), & 2926 'lat', NF90_DOUBLE, id_var_lat_xy(av,0), & 2927 'degrees_north', 'latitude', 000, 000, 000 ) 2928 CALL netcdf_create_var( id_set_xy(av), & 2929 (/ id_dim_xu_xy(av), id_dim_y_xy(av) /), & 2930 'lonu', NF90_DOUBLE, id_var_lon_xy(av,1), & 2931 'degrees_east', 'longitude', 000, 000, 000 ) 2932 CALL netcdf_create_var( id_set_xy(av), & 2933 (/ id_dim_xu_xy(av), id_dim_y_xy(av) /), & 2934 'latu', NF90_DOUBLE, id_var_lat_xy(av,1), & 2935 'degrees_north', 'latitude', 000, 000, 000 ) 2936 CALL netcdf_create_var( id_set_xy(av), & 2937 (/ id_dim_x_xy(av), id_dim_yv_xy(av) /), & 2938 'lonv', NF90_DOUBLE, id_var_lon_xy(av,2), & 2939 'degrees_east', 'longitude', 000, 000, 000 ) 2940 CALL netcdf_create_var( id_set_xy(av), & 2941 (/ id_dim_x_xy(av), id_dim_yv_xy(av) /), & 2942 'latv', NF90_DOUBLE, id_var_lat_xy(av,2), & 2943 'degrees_north', 'latitude', 000, 000, 000 ) 2944 ! 2945 !-- Define coordinate-reference system 2946 CALL netcdf_create_crs( id_set_xy(av), 000 ) 2722 2947 ! 2723 2948 !-- In case of non-flat topography define 2d-arrays containing the height … … 3199 3424 ENDIF 3200 3425 ! 3426 !-- Write lon and lat data. Only for parallel output. 3427 IF ( netcdf_data_format > 4 ) THEN 3428 3429 ALLOCATE( lat(nxl:nxr,nys:nyn) ) 3430 ALLOCATE( lon(nxl:nxr,nys:nyn) ) 3431 cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp ) 3432 sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp ) 3433 3434 DO k = 0, 2 3435 ! 3436 !-- Scalar grid points 3437 IF ( k == 0 ) THEN 3438 shift_x = 0.5 ; shift_y = 0.5 3439 ! 3440 !-- u grid points 3441 ELSEIF ( k == 1 ) THEN 3442 shift_x = 0.0 ; shift_y = 0.5 3443 ! 3444 !-- v grid points 3445 ELSEIF ( k == 2 ) THEN 3446 shift_x = 0.5 ; shift_y = 0.0 3447 ENDIF 3448 3449 DO j = nys, nyn 3450 DO i = nxl, nxr 3451 eutm = init_model%origin_x & 3452 + cos_ra * ( i + shift_x ) * dx & 3453 + sin_ra * ( j + shift_y ) * dy 3454 nutm = init_model%origin_y & 3455 - sin_ra * ( i + shift_x ) * dx & 3456 + cos_ra * ( j + shift_y ) * dy 3457 3458 CALL convert_utm_to_geographic( crs_list, & 3459 eutm, nutm, & 3460 lon(i,j), lat(i,j) ) 3461 ENDDO 3462 ENDDO 3463 3464 nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_lon_xy(av,k), & 3465 lon, start = (/ nxl+1, nys+1 /), & 3466 count = (/ nxr-nxl+1, nyn-nys+1 /) ) 3467 CALL netcdf_handle_error( 'netcdf_define_header', 556 ) 3468 3469 nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_lat_xy(av,k), & 3470 lat, start = (/ nxl+1, nys+1 /), & 3471 count = (/ nxr-nxl+1, nyn-nys+1 /) ) 3472 CALL netcdf_handle_error( 'netcdf_define_header', 556 ) 3473 ENDDO 3474 3475 DEALLOCATE( lat ) 3476 DEALLOCATE( lon ) 3477 3478 ENDIF 3479 ! 3201 3480 !-- In case of non-flat topography write height information. Only for 3202 3481 !-- parallel netcdf output. … … 3595 3874 (/ id_dim_x_xz(av) /), & 3596 3875 'E_UTM', NF90_DOUBLE, id_var_eutm_xz(av,0), & 3597 'm', ' ', 000, 000, 000 )3876 'm', 'projection_x_coordinate', 000, 000, 000 ) 3598 3877 CALL netcdf_create_var( id_set_xz(av), & 3599 3878 (/ id_dim_y_xz(av) /), & 3600 3879 'N_UTM', NF90_DOUBLE, id_var_nutm_xz(av,0), & 3601 'm', ' ', 000, 000, 000 )3880 'm', 'projection_y_coordinate', 000, 000, 000 ) 3602 3881 CALL netcdf_create_var( id_set_xz(av), & 3603 3882 (/ id_dim_xu_xz(av) /), & 3604 3883 'Eu_UTM', NF90_DOUBLE, id_var_eutm_xz(av,1), & 3605 'm', ' ', 000, 000, 000 )3884 'm', 'projection_x_coordinate', 000, 000, 000 ) 3606 3885 CALL netcdf_create_var( id_set_xz(av), & 3607 3886 (/ id_dim_y_xz(av) /), & 3608 3887 'Nu_UTM', NF90_DOUBLE, id_var_nutm_xz(av,1), & 3609 'm', ' ', 000, 000, 000 )3888 'm', 'projection_y_coordinate', 000, 000, 000 ) 3610 3889 CALL netcdf_create_var( id_set_xz(av), & 3611 3890 (/ id_dim_x_xz(av) /), & 3612 3891 'Ev_UTM', NF90_DOUBLE, id_var_eutm_xz(av,2), & 3613 'm', ' ', 000, 000, 000 )3892 'm', 'projection_x_coordinate', 000, 000, 000 ) 3614 3893 CALL netcdf_create_var( id_set_xz(av), & 3615 3894 (/ id_dim_yv_xz(av) /), & 3616 3895 'Nv_UTM', NF90_DOUBLE, id_var_nutm_xz(av,2), & 3617 'm', ' ', 000, 000, 000 )3896 'm', 'projection_y_coordinate', 000, 000, 000 ) 3618 3897 ELSE 3619 3898 CALL netcdf_create_var( id_set_xz(av), & 3620 3899 (/ id_dim_x_xz(av), id_dim_y_xz(av) /), & 3621 3900 'E_UTM', NF90_DOUBLE, id_var_eutm_xz(av,0), & 3622 'm', ' ', 000, 000, 000 )3901 'm', 'projection_x_coordinate', 000, 000, 000 ) 3623 3902 CALL netcdf_create_var( id_set_xz(av), & 3624 3903 (/ id_dim_x_xz(av), id_dim_y_xz(av) /), & 3625 3904 'N_UTM', NF90_DOUBLE, id_var_nutm_xz(av,0), & 3626 'm', ' ', 000, 000, 000 )3905 'm', 'projection_y_coordinate', 000, 000, 000 ) 3627 3906 CALL netcdf_create_var( id_set_xz(av), & 3628 3907 (/ id_dim_xu_xz(av), id_dim_y_xz(av) /), & 3629 3908 'Eu_UTM', NF90_DOUBLE, id_var_eutm_xz(av,1), & 3630 'm', ' ', 000, 000, 000 )3909 'm', 'projection_x_coordinate', 000, 000, 000 ) 3631 3910 CALL netcdf_create_var( id_set_xz(av), & 3632 3911 (/ id_dim_xu_xz(av), id_dim_y_xz(av) /), & 3633 3912 'Nu_UTM', NF90_DOUBLE, id_var_nutm_xz(av,1), & 3634 'm', ' ', 000, 000, 000 )3913 'm', 'projection_y_coordinate', 000, 000, 000 ) 3635 3914 CALL netcdf_create_var( id_set_xz(av), & 3636 3915 (/ id_dim_x_xz(av), id_dim_yv_xz(av) /), & 3637 3916 'Ev_UTM', NF90_DOUBLE, id_var_eutm_xz(av,2), & 3638 'm', ' ', 000, 000, 000 )3917 'm', 'projection_x_coordinate', 000, 000, 000 ) 3639 3918 CALL netcdf_create_var( id_set_xz(av), & 3640 3919 (/ id_dim_x_xz(av), id_dim_yv_xz(av) /), & 3641 3920 'Nv_UTM', NF90_DOUBLE, id_var_nutm_xz(av,2), & 3642 'm', '', 000, 000, 000 ) 3643 ENDIF 3921 'm', 'projection_y_coordinate', 000, 000, 000 ) 3922 ENDIF 3923 ! 3924 !-- Define geographic coordinates 3925 CALL netcdf_create_var( id_set_xz(av), & 3926 (/ id_dim_x_xz(av), id_dim_y_xz(av) /), & 3927 'lon', NF90_DOUBLE, id_var_lon_xz(av,0), & 3928 'degrees_east', 'longitude', 000, 000, 000 ) 3929 CALL netcdf_create_var( id_set_xz(av), & 3930 (/ id_dim_x_xz(av), id_dim_y_xz(av) /), & 3931 'lat', NF90_DOUBLE, id_var_lat_xz(av,0), & 3932 'degrees_north', 'latitude', 000, 000, 000 ) 3933 CALL netcdf_create_var( id_set_xz(av), & 3934 (/ id_dim_xu_xz(av), id_dim_y_xz(av) /), & 3935 'lonu', NF90_DOUBLE, id_var_lon_xz(av,1), & 3936 'degrees_east', 'longitude', 000, 000, 000 ) 3937 CALL netcdf_create_var( id_set_xz(av), & 3938 (/ id_dim_xu_xz(av), id_dim_y_xz(av) /), & 3939 'latu', NF90_DOUBLE, id_var_lat_xz(av,1), & 3940 'degrees_north', 'latitude', 000, 000, 000 ) 3941 CALL netcdf_create_var( id_set_xz(av), & 3942 (/ id_dim_x_xz(av), id_dim_yv_xz(av) /), & 3943 'lonv', NF90_DOUBLE, id_var_lon_xz(av,2), & 3944 'degrees_east', 'longitude', 000, 000, 000 ) 3945 CALL netcdf_create_var( id_set_xz(av), & 3946 (/ id_dim_x_xz(av), id_dim_yv_xz(av) /), & 3947 'latv', NF90_DOUBLE, id_var_lat_xz(av,2), & 3948 'degrees_north', 'latitude', 000, 000, 000 ) 3949 ! 3950 !-- Define coordinate-reference system 3951 CALL netcdf_create_crs( id_set_xz(av), 000 ) 3644 3952 3645 3953 IF ( land_surface ) THEN … … 4080 4388 DEALLOCATE( netcdf_data_2d ) 4081 4389 ENDIF 4390 ! 4391 !-- Write lon and lat data 4392 ALLOCATE( lat(0:nx,1:ns) ) 4393 ALLOCATE( lon(0:nx,1:ns) ) 4394 cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp ) 4395 sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp ) 4396 4397 DO k = 0, 2 4398 ! 4399 !-- Scalar grid points 4400 IF ( k == 0 ) THEN 4401 shift_x = 0.5 ; shift_y = 0.5 4402 ! 4403 !-- u grid points 4404 ELSEIF ( k == 1 ) THEN 4405 shift_x = 0.0 ; shift_y = 0.5 4406 ! 4407 !-- v grid points 4408 ELSEIF ( k == 2 ) THEN 4409 shift_x = 0.5 ; shift_y = 0.0 4410 ENDIF 4411 4412 DO j = 1, ns 4413 IF( section(j,2) == -1 ) THEN 4414 lat(:,j) = -90.0_wp ! section averaged along y 4415 lon(:,j) = -180.0_wp ! section averaged along y 4416 ELSE 4417 DO i = 0, nx 4418 eutm = init_model%origin_x & 4419 + cos_ra * ( i + shift_x ) * dx & 4420 + sin_ra * ( section(j,2) + shift_y ) * dy 4421 nutm = init_model%origin_y & 4422 - sin_ra * ( i + shift_x ) * dx & 4423 + cos_ra * ( section(j,2) + shift_y ) * dy 4424 4425 CALL convert_utm_to_geographic( crs_list, & 4426 eutm, nutm, & 4427 lon(i,j), lat(i,j) ) 4428 ENDDO 4429 ENDIF 4430 ENDDO 4431 4432 nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_lon_xz(av,k), & 4433 lon, start = (/ 1, 1 /), & 4434 count = (/ nx+1, ns /) ) 4435 CALL netcdf_handle_error( 'netcdf_define_header', 556 ) 4436 4437 nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_lat_xz(av,k), & 4438 lat, start = (/ 1, 1 /), & 4439 count = (/ nx+1, ns /) ) 4440 CALL netcdf_handle_error( 'netcdf_define_header', 556 ) 4441 ENDDO 4442 4443 DEALLOCATE( lat ) 4444 DEALLOCATE( lon ) 4082 4445 4083 4446 ENDIF … … 4444 4807 (/ id_dim_x_yz(av) /), & 4445 4808 'E_UTM', NF90_DOUBLE, id_var_eutm_yz(av,0), & 4446 'm', ' ', 000, 000, 000 )4809 'm', 'projection_x_coordinate', 000, 000, 000 ) 4447 4810 CALL netcdf_create_var( id_set_yz(av), & 4448 4811 (/ id_dim_y_yz(av) /), & 4449 4812 'N_UTM', NF90_DOUBLE, id_var_nutm_yz(av,0), & 4450 'm', ' ', 000, 000, 000 )4813 'm', 'projection_y_coordinate', 000, 000, 000 ) 4451 4814 CALL netcdf_create_var( id_set_yz(av), & 4452 4815 (/ id_dim_xu_yz(av) /), & 4453 4816 'Eu_UTM', NF90_DOUBLE, id_var_eutm_yz(av,1), & 4454 'm', ' ', 000, 000, 000 )4817 'm', 'projection_x_coordinate', 000, 000, 000 ) 4455 4818 CALL netcdf_create_var( id_set_yz(av), & 4456 4819 (/ id_dim_y_yz(av) /), & 4457 4820 'Nu_UTM', NF90_DOUBLE, id_var_nutm_yz(av,1), & 4458 'm', ' ', 000, 000, 000 )4821 'm', 'projection_y_coordinate', 000, 000, 000 ) 4459 4822 CALL netcdf_create_var( id_set_yz(av), & 4460 4823 (/ id_dim_x_yz(av) /), & 4461 4824 'Ev_UTM', NF90_DOUBLE, id_var_eutm_yz(av,2), & 4462 'm', ' ', 000, 000, 000 )4825 'm', 'projection_x_coordinate', 000, 000, 000 ) 4463 4826 CALL netcdf_create_var( id_set_yz(av), & 4464 4827 (/ id_dim_yv_yz(av) /), & 4465 4828 'Nv_UTM', NF90_DOUBLE, id_var_nutm_yz(av,2), & 4466 'm', ' ', 000, 000, 000 )4829 'm', 'projection_y_coordinate', 000, 000, 000 ) 4467 4830 ELSE 4468 4831 CALL netcdf_create_var( id_set_yz(av), & 4469 4832 (/ id_dim_x_yz(av), id_dim_y_yz(av) /), & 4470 4833 'E_UTM', NF90_DOUBLE, id_var_eutm_yz(av,0), & 4471 'm', ' ', 000, 000, 000 )4834 'm', 'projection_x_coordinate', 000, 000, 000 ) 4472 4835 CALL netcdf_create_var( id_set_yz(av), & 4473 4836 (/ id_dim_x_yz(av), id_dim_y_yz(av) /), & 4474 4837 'N_UTM', NF90_DOUBLE, id_var_nutm_yz(av,0), & 4475 'm', ' ', 000, 000, 000 )4838 'm', 'projection_y_coordinate', 000, 000, 000 ) 4476 4839 CALL netcdf_create_var( id_set_yz(av), & 4477 4840 (/ id_dim_xu_yz(av), id_dim_y_yz(av) /), & 4478 4841 'Eu_UTM', NF90_DOUBLE, id_var_eutm_yz(av,1), & 4479 'm', ' ', 000, 000, 000 )4842 'm', 'projection_x_coordinate', 000, 000, 000 ) 4480 4843 CALL netcdf_create_var( id_set_yz(av), & 4481 4844 (/ id_dim_xu_yz(av), id_dim_y_yz(av) /), & 4482 4845 'Nu_UTM', NF90_DOUBLE, id_var_nutm_yz(av,1), & 4483 'm', ' ', 000, 000, 000 )4846 'm', 'projection_y_coordinate', 000, 000, 000 ) 4484 4847 CALL netcdf_create_var( id_set_yz(av), & 4485 4848 (/ id_dim_x_yz(av), id_dim_yv_yz(av) /), & 4486 4849 'Ev_UTM', NF90_DOUBLE, id_var_eutm_yz(av,2), & 4487 'm', ' ', 000, 000, 000 )4850 'm', 'projection_x_coordinate', 000, 000, 000 ) 4488 4851 CALL netcdf_create_var( id_set_yz(av), & 4489 4852 (/ id_dim_x_yz(av), id_dim_yv_yz(av) /), & 4490 4853 'Nv_UTM', NF90_DOUBLE, id_var_nutm_yz(av,2), & 4491 'm', '', 000, 000, 000 ) 4492 ENDIF 4854 'm', 'projection_y_coordinate', 000, 000, 000 ) 4855 ENDIF 4856 ! 4857 !-- Define geographic coordinates 4858 CALL netcdf_create_var( id_set_yz(av), & 4859 (/ id_dim_x_yz(av), id_dim_y_yz(av) /), & 4860 'lon', NF90_DOUBLE, id_var_lon_yz(av,0), & 4861 'degrees_east', 'longitude', 000, 000, 000 ) 4862 CALL netcdf_create_var( id_set_yz(av), & 4863 (/ id_dim_x_yz(av), id_dim_y_yz(av) /), & 4864 'lat', NF90_DOUBLE, id_var_lat_yz(av,0), & 4865 'degrees_north', 'latitude', 000, 000, 000 ) 4866 CALL netcdf_create_var( id_set_yz(av), & 4867 (/ id_dim_xu_yz(av), id_dim_y_yz(av) /), & 4868 'lonu', NF90_DOUBLE, id_var_lon_yz(av,1), & 4869 'degrees_east', 'longitude', 000, 000, 000 ) 4870 CALL netcdf_create_var( id_set_yz(av), & 4871 (/ id_dim_xu_yz(av), id_dim_y_yz(av) /), & 4872 'latu', NF90_DOUBLE, id_var_lat_yz(av,1), & 4873 'degrees_north', 'latitude', 000, 000, 000 ) 4874 CALL netcdf_create_var( id_set_yz(av), & 4875 (/ id_dim_x_yz(av), id_dim_yv_yz(av) /), & 4876 'lonv', NF90_DOUBLE, id_var_lon_yz(av,2), & 4877 'degrees_east', 'longitude', 000, 000, 000 ) 4878 CALL netcdf_create_var( id_set_yz(av), & 4879 (/ id_dim_x_yz(av), id_dim_yv_yz(av) /), & 4880 'latv', NF90_DOUBLE, id_var_lat_yz(av,2), & 4881 'degrees_north', 'latitude', 000, 000, 000 ) 4882 ! 4883 !-- Define coordinate-reference system 4884 CALL netcdf_create_crs( id_set_yz(av), 000 ) 4493 4885 4494 4886 IF ( land_surface ) THEN … … 4916 5308 DEALLOCATE( netcdf_data_2d ) 4917 5309 ENDIF 5310 ! 5311 !-- Write lon and lat data 5312 ALLOCATE( lat(1:ns,0:ny) ) 5313 ALLOCATE( lon(1:ns,0:ny) ) 5314 cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp ) 5315 sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp ) 5316 5317 DO k = 0, 2 5318 ! 5319 !-- Scalar grid points 5320 IF ( k == 0 ) THEN 5321 shift_x = 0.5 ; shift_y = 0.5 5322 ! 5323 !-- u grid points 5324 ELSEIF ( k == 1 ) THEN 5325 shift_x = 0.0 ; shift_y = 0.5 5326 ! 5327 !-- v grid points 5328 ELSEIF ( k == 2 ) THEN 5329 shift_x = 0.5 ; shift_y = 0.0 5330 ENDIF 5331 5332 DO j = 0, ny 5333 DO i = 1, ns 5334 IF( section(i,3) == -1 ) THEN 5335 lat(i,:) = -90.0_wp ! section averaged along x 5336 lon(i,:) = -180.0_wp ! section averaged along x 5337 ELSE 5338 eutm = init_model%origin_x & 5339 + cos_ra * ( section(i,3) + shift_x ) * dx & 5340 + sin_ra * ( j + shift_y ) * dy 5341 nutm = init_model%origin_y & 5342 - sin_ra * ( section(i,3) + shift_x ) * dx & 5343 + cos_ra * ( j + shift_y ) * dy 5344 5345 CALL convert_utm_to_geographic( crs_list, & 5346 eutm, nutm, & 5347 lon(i,j), lat(i,j) ) 5348 ENDIF 5349 ENDDO 5350 ENDDO 5351 5352 nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_lon_yz(av,k), & 5353 lon, start = (/ 1, 1 /), & 5354 count = (/ ns, ny+1 /) ) 5355 CALL netcdf_handle_error( 'netcdf_define_header', 556 ) 5356 5357 nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_lat_yz(av,k), & 5358 lat, start = (/ 1, 1 /), & 5359 count = (/ ns, ny+1 /) ) 5360 CALL netcdf_handle_error( 'netcdf_define_header', 556 ) 5361 ENDDO 5362 5363 DEALLOCATE( lat ) 5364 DEALLOCATE( lon ) 4918 5365 4919 5366 ENDIF … … 6896 7343 END SUBROUTINE netcdf_create_var 6897 7344 7345 7346 !------------------------------------------------------------------------------! 7347 ! Description: 7348 ! ------------ 7349 !> Create a variable holding the coordinate-reference-system information. 7350 !------------------------------------------------------------------------------! 7351 SUBROUTINE netcdf_create_crs( ncid, error_no ) 7352 7353 #if defined( __netcdf ) 7354 IMPLICIT NONE 7355 7356 INTEGER, INTENT(IN) :: error_no !< error number 7357 INTEGER, INTENT(IN) :: ncid !< file id 7358 INTEGER :: var_id !< variable id 7359 7360 ! 7361 !-- Define variable 7362 nc_stat = NF90_DEF_VAR( ncid, 'crs', NF90_INT, VARID = var_id ) 7363 CALL netcdf_handle_error( 'netcdf_create_crs', error_no ) 7364 ! 7365 !-- Set attributes 7366 nc_stat = NF90_PUT_ATT( ncid, var_id, 'epsg_code', & 7367 coord_ref_sys%epsg_code ) 7368 CALL netcdf_handle_error( 'netcdf_create_crs', error_no ) 7369 7370 nc_stat = NF90_PUT_ATT( ncid, var_id, 'false_easting', & 7371 coord_ref_sys%false_easting ) 7372 CALL netcdf_handle_error( 'netcdf_create_crs', error_no ) 7373 7374 nc_stat = NF90_PUT_ATT( ncid, var_id, 'false_northing', & 7375 coord_ref_sys%false_northing ) 7376 CALL netcdf_handle_error( 'netcdf_create_crs', error_no ) 7377 7378 nc_stat = NF90_PUT_ATT( ncid, var_id, 'grid_mapping_name', & 7379 coord_ref_sys%grid_mapping_name ) 7380 CALL netcdf_handle_error( 'netcdf_create_crs', error_no ) 7381 7382 nc_stat = NF90_PUT_ATT( ncid, var_id, 'inverse_flattening', & 7383 coord_ref_sys%inverse_flattening ) 7384 CALL netcdf_handle_error( 'netcdf_create_crs', error_no ) 7385 7386 nc_stat = NF90_PUT_ATT( ncid, var_id, 'latitude_of_projection_origin', & 7387 coord_ref_sys%latitude_of_projection_origin ) 7388 CALL netcdf_handle_error( 'netcdf_create_crs', error_no ) 7389 7390 nc_stat = NF90_PUT_ATT( ncid, var_id, 'long_name', & 7391 coord_ref_sys%long_name ) 7392 CALL netcdf_handle_error( 'netcdf_create_crs', error_no ) 7393 7394 nc_stat = NF90_PUT_ATT( ncid, var_id, 'longitude_of_central_meridian', & 7395 coord_ref_sys%longitude_of_central_meridian ) 7396 CALL netcdf_handle_error( 'netcdf_create_crs', error_no ) 7397 7398 nc_stat = NF90_PUT_ATT( ncid, var_id, 'longitude_of_prime_meridian', & 7399 coord_ref_sys%longitude_of_prime_meridian ) 7400 CALL netcdf_handle_error( 'netcdf_create_crs', error_no ) 7401 7402 nc_stat = NF90_PUT_ATT( ncid, var_id, 'scale_factor_at_central_meridian', & 7403 coord_ref_sys%scale_factor_at_central_meridian ) 7404 CALL netcdf_handle_error( 'netcdf_create_crs', error_no ) 7405 7406 nc_stat = NF90_PUT_ATT( ncid, var_id, 'semi_major_axis', & 7407 coord_ref_sys%semi_major_axis ) 7408 CALL netcdf_handle_error( 'netcdf_create_crs', error_no ) 7409 7410 nc_stat = NF90_PUT_ATT( ncid, var_id, 'units', & 7411 coord_ref_sys%units ) 7412 CALL netcdf_handle_error( 'netcdf_create_crs', error_no ) 7413 7414 #endif 7415 END SUBROUTINE netcdf_create_crs 7416 7417 7418 !------------------------------------------------------------------------------! 7419 ! Description: 7420 ! ------------ 7421 !> Convert UTM coordinates into geographic latitude and longitude. Conversion 7422 !> is based on the work of KrÃŒger (1912) DOI: 10.2312/GFZ.b103-krueger28 7423 !> and Karney (2013) DOI: 10.1007/s00190-012-0578-z 7424 !> Based on a JavaScript of the geodesy function library written by chrisveness 7425 !> https://github.com/chrisveness/geodesy 7426 !------------------------------------------------------------------------------! 7427 SUBROUTINE convert_utm_to_geographic( crs, eutm, nutm, lon, lat ) 7428 7429 USE basic_constants_and_equations_mod, & 7430 ONLY: pi 7431 7432 IMPLICIT NONE 7433 7434 REAL(wp), INTENT(in) :: eutm !< easting (UTM) 7435 REAL(wp), INTENT(out) :: lat !< geographic latitude in degree 7436 REAL(wp), INTENT(out) :: lon !< geographic longitude in degree 7437 REAL(wp), INTENT(in) :: nutm !< northing (UTM) 7438 7439 REAL(wp) :: a !< 2*pi*a is the circumference of a meridian 7440 REAL(wp) :: cos_eta_s !< cos(eta_s) 7441 REAL(wp) :: delta_i !< 7442 REAL(wp) :: delta_tau_i !< 7443 REAL(wp) :: e !< eccentricity 7444 REAL(wp) :: eta !< 7445 REAL(wp) :: eta_s !< 7446 REAL(wp) :: j !< loop index 7447 REAL(wp) :: n !< 3rd flattening 7448 REAL(wp) :: n2 !< n^2 7449 REAL(wp) :: n3 !< n^3 7450 REAL(wp) :: n4 !< n^4 7451 REAL(wp) :: n5 !< n^5 7452 REAL(wp) :: n6 !< n^6 7453 REAL(wp) :: nu !< 7454 REAL(wp) :: nu_s !< 7455 REAL(wp) :: sin_eta_s !< sin(eta_s) 7456 REAL(wp) :: sinh_nu_s !< sinush(nu_s) 7457 REAL(wp) :: tau_i !< 7458 REAL(wp) :: tau_i_s !< 7459 REAL(wp) :: tau_s !< 7460 REAL(wp) :: x !< adjusted easting 7461 REAL(wp) :: y !< adjusted northing 7462 7463 REAL(wp), DIMENSION(6) :: beta !< 6th order KrÃŒger expressions 7464 7465 REAL(wp), DIMENSION(8), INTENT(in) :: crs !< coordinate reference system, consists of 7466 !< (/semi_major_axis, 7467 !< inverse_flattening, 7468 !< longitude_of_prime_meridian, 7469 !< longitude_of_central_meridian, 7470 !< scale_factor_at_central_meridian, 7471 !< latitude_of_projection_origin, 7472 !< false_easting, 7473 !< false_northing /) 7474 7475 x = eutm - crs(7) ! remove false easting 7476 y = nutm - crs(8) ! remove false northing 7477 7478 !-- from Karney 2011 Eq 15-22, 36: 7479 e = SQRT( 1.0_wp / crs(2) * ( 2.0_wp - 1.0_wp / crs(2) ) ) 7480 n = 1.0_wp / crs(2) / ( 2.0_wp - 1.0_wp / crs(2) ) 7481 n2 = n * n 7482 n3 = n * n2 7483 n4 = n * n3 7484 n5 = n * n4 7485 n6 = n * n5 7486 7487 a = crs(1) / ( 1.0_wp + n ) * ( 1.0_wp + 0.25_wp * n2 & 7488 + 0.015625_wp * n4 & 7489 + 3.90625E-3_wp * n6 ) 7490 7491 nu = x / ( crs(5) * a ) 7492 eta = y / ( crs(5) * a ) 7493 7494 !-- According to KrÃŒger (1912), eq. 26* 7495 beta(1) = 0.5_wp * n & 7496 - 2.0_wp / 3.0_wp * n2 & 7497 + 37.0_wp / 96.0_wp * n3 & 7498 - 1.0_wp / 360.0_wp * n4 & 7499 - 81.0_wp / 512.0_wp * n5 & 7500 + 96199.0_wp / 604800.0_wp * n6 7501 7502 beta(2) = 1.0_wp / 48.0_wp * n2 & 7503 + 1.0_wp / 15.0_wp * n3 & 7504 - 437.0_wp / 1440.0_wp * n4 & 7505 + 46.0_wp / 105.0_wp * n5 & 7506 - 1118711.0_wp / 3870720.0_wp * n6 7507 7508 beta(3) = 17.0_wp / 480.0_wp * n3 & 7509 - 37.0_wp / 840.0_wp * n4 & 7510 - 209.0_wp / 4480.0_wp * n5 & 7511 + 5569.0_wp / 90720.0_wp * n6 7512 7513 beta(4) = 4397.0_wp / 161280.0_wp * n4 & 7514 - 11.0_wp / 504.0_wp * n5 & 7515 - 830251.0_wp / 7257600.0_wp * n6 7516 7517 beta(5) = 4583.0_wp / 161280.0_wp * n5 & 7518 - 108847.0_wp / 3991680.0_wp * n6 7519 7520 beta(6) = 20648693.0_wp / 638668800.0_wp * n6 7521 7522 eta_s = eta 7523 nu_s = nu 7524 DO j = 1, 6 7525 eta_s = eta_s - beta(j) * SIN(2.0_wp * j * eta) * COSH(2.0_wp * j * nu) 7526 nu_s = nu_s - beta(j) * COS(2.0_wp * j * eta) * SINH(2.0_wp * j * nu) 7527 ENDDO 7528 7529 sinh_nu_s = SINH( nu_s ) 7530 sin_eta_s = SIN( eta_s ) 7531 cos_eta_s = COS( eta_s ) 7532 7533 tau_s = sin_eta_s / SQRT( sinh_nu_s**2 + cos_eta_s**2 ) 7534 7535 tau_i = tau_s 7536 delta_tau_i = 1.0_wp 7537 7538 DO WHILE ( ABS( delta_tau_i ) > 1.0E-12_wp ) 7539 7540 delta_i = SINH( e * ATANH( e * tau_i / SQRT( 1.0_wp + tau_i**2 ) ) ) 7541 7542 tau_i_s = tau_i * SQRT( 1.0_wp + delta_i**2 ) & 7543 - delta_i * SQRT( 1.0_wp + tau_i**2 ) 7544 7545 delta_tau_i = ( tau_s - tau_i_s ) / SQRT( 1.0_wp + tau_i_s**2 ) & 7546 * ( 1.0_wp + ( 1.0_wp - e**2 ) * tau_i**2 ) & 7547 / ( ( 1.0_wp - e**2 ) * SQRT( 1.0_wp + tau_i**2 ) ) 7548 7549 tau_i = tau_i + delta_tau_i 7550 7551 ENDDO 7552 7553 lat = ATAN( tau_i ) / pi * 180.0_wp 7554 lon = ATAN2( sinh_nu_s, cos_eta_s ) / pi * 180.0_wp + crs(4) 7555 7556 END SUBROUTINE convert_utm_to_geographic 7557 6898 7558 END MODULE netcdf_interface
Note: See TracChangeset
for help on using the changeset viewer.
