Ignore:
Timestamp:
Oct 30, 2018 3:04:11 PM (3 years ago)
Author:
gronemeier
Message:

add longitude and latitude to output files; update run_control files for tests

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/netcdf_interface_mod.f90

    r3448 r3459  
    2020! Current revisions:
    2121! ------------------
    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.
    2326!
    2427! Former revisions:
     
    311314!> @todo timeseries and profile output still needs to be rewritten to allow
    312315!>       modularization
     316!> @todo output 2d UTM coordinates without global arrays
     317!> @todo output longitude/latitude also with non-parallel output (3d and xy)
    313318!------------------------------------------------------------------------------!
    314319 MODULE netcdf_interface
     
    324329    USE mas_global_attributes,                                                 &
    325330        ONLY:  dim_size_agtnum
     331
     332    USE netcdf_data_input_mod,                                                 &
     333        ONLY: coord_ref_sys
    326334
    327335    PRIVATE
     
    472480                                         id_var_eutm_xy, id_var_nutm_xy, &
    473481                                         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
    476488
    477489    INTEGER ::  netcdf_data_format = 2  !< NetCDF3 64bit offset format
     
    520532    INTEGER(iwp), DIMENSION(1:max_masks,0:1,0:2) ::  id_var_eutm_mask, &
    521533                                                     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
    522537
    523538    INTEGER(iwp), DIMENSION(1:max_masks,0:1,100) ::  id_var_domask
     
    624639
    625640    USE netcdf_data_input_mod,                                                 &
    626         ONLY: init_model
     641        ONLY:  init_model
    627642
    628643    USE ocean_mod,                                                             &
     
    720735
    721736    REAL(wp) ::  cos_ra                                      !< cosine of rotation_angle
     737    REAL(wp) ::  eutm                                        !< easting (UTM)
     738    REAL(wp) ::  nutm                                        !< northing (UTM)
    722739    REAL(wp) ::  shift_x                                     !< shift of x coordinate
    723740    REAL(wp) ::  shift_y                                     !< shift of y coordinate
     
    725742
    726743    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
    727745
    728746    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  netcdf_data    !<
    729747    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  netcdf_data_2d !<
     748    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lat            !< latitude
     749    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lon            !< longitude
    730750
    731751
     
    822842
    823843!
     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!
    824855!-- Determine the mode to be processed
    825856    IF ( extend )  THEN
     
    962993                                 (/ id_dim_x_mask(mid,av) /),      &
    963994                                 'E_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,0),  &
    964                                  'm', '', 000, 000, 000 )
     995                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    965996             CALL netcdf_create_var( id_set_mask(mid,av), &
    966997                                 (/ id_dim_y_mask(mid,av) /),      &
    967998                                 'N_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,0),  &
    968                                  'm', '', 000, 000, 000 )
     999                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    9691000             CALL netcdf_create_var( id_set_mask(mid,av), &
    9701001                                 (/ id_dim_xu_mask(mid,av) /),     &
    9711002                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,1), &
    972                                  'm', '', 000, 000, 000 )
     1003                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    9731004             CALL netcdf_create_var( id_set_mask(mid,av), &
    9741005                                 (/ id_dim_y_mask(mid,av) /),     &
    9751006                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,1), &
    976                                  'm', '', 000, 000, 000 )
     1007                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    9771008             CALL netcdf_create_var( id_set_mask(mid,av), &
    9781009                                 (/ id_dim_x_mask(mid,av) /),     &
    9791010                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,2), &
    980                                  'm', '', 000, 000, 000 )
     1011                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    9811012             CALL netcdf_create_var( id_set_mask(mid,av), &
    9821013                                 (/ id_dim_yv_mask(mid,av) /),     &
    9831014                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,2), &
    984                                  'm', '', 000, 000, 000 )
     1015                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    9851016          ELSE
    9861017             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 )
    9901021             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 )
    9941025             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 )
    9981029             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 )
    10021033             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 )
    10061037             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 )
    10111071!
    10121072!--       In case of non-flat topography define 2d-arrays containing the height
     
    13641424             DEALLOCATE( netcdf_data_2d )
    13651425          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 )
    13671480!
    13681481!--       Write zu and zw data (vertical axes)
     
    17031816                                 (/ id_dim_x_3d(av) /),      &
    17041817                                 'E_UTM', NF90_DOUBLE, id_var_eutm_3d(av,0),  &
    1705                                  'm', '', 000, 000, 000 )
     1818                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    17061819             CALL netcdf_create_var( id_set_3d(av), &
    17071820                                 (/ id_dim_y_3d(av) /),      &
    17081821                                 'N_UTM', NF90_DOUBLE, id_var_nutm_3d(av,0),  &
    1709                                  'm', '', 000, 000, 000 )
     1822                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    17101823             CALL netcdf_create_var( id_set_3d(av), &
    17111824                                 (/ id_dim_xu_3d(av) /),     &
    17121825                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_3d(av,1), &
    1713                                  'm', '', 000, 000, 000 )
     1826                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    17141827             CALL netcdf_create_var( id_set_3d(av), &
    17151828                                 (/ id_dim_y_3d(av) /),     &
    17161829                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_3d(av,1), &
    1717                                  'm', '', 000, 000, 000 )
     1830                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    17181831             CALL netcdf_create_var( id_set_3d(av), &
    17191832                                 (/ id_dim_x_3d(av) /),     &
    17201833                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_3d(av,2), &
    1721                                  'm', '', 000, 000, 000 )
     1834                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    17221835             CALL netcdf_create_var( id_set_3d(av), &
    17231836                                 (/ id_dim_yv_3d(av) /),     &
    17241837                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_3d(av,2), &
    1725                                  'm', '', 000, 000, 000 )
     1838                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    17261839          ELSE
    17271840             CALL netcdf_create_var( id_set_3d(av), &
    17281841                                 (/ id_dim_x_3d(av), id_dim_y_3d(av) /),      &
    17291842                                 'E_UTM', NF90_DOUBLE, id_var_eutm_3d(av,0),  &
    1730                                  'm', '', 000, 000, 000 )
     1843                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    17311844             CALL netcdf_create_var( id_set_3d(av), &
    17321845                                 (/ id_dim_x_3d(av), id_dim_y_3d(av) /),      &
    17331846                                 'N_UTM', NF90_DOUBLE, id_var_nutm_3d(av,0),  &
    1734                                  'm', '', 000, 000, 000 )
     1847                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    17351848             CALL netcdf_create_var( id_set_3d(av), &
    17361849                                 (/ id_dim_xu_3d(av), id_dim_y_3d(av) /),     &
    17371850                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_3d(av,1), &
    1738                                  'm', '', 000, 000, 000 )
     1851                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    17391852             CALL netcdf_create_var( id_set_3d(av), &
    17401853                                 (/ id_dim_xu_3d(av), id_dim_y_3d(av) /),     &
    17411854                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_3d(av,1), &
    1742                                  'm', '', 000, 000, 000 )
     1855                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    17431856             CALL netcdf_create_var( id_set_3d(av), &
    17441857                                 (/ id_dim_x_3d(av), id_dim_yv_3d(av) /),     &
    17451858                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_3d(av,2), &
    1746                                  'm', '', 000, 000, 000 )
     1859                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    17471860             CALL netcdf_create_var( id_set_3d(av), &
    17481861                                 (/ id_dim_x_3d(av), id_dim_yv_3d(av) /),     &
    17491862                                 '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 )
    17521894!
    17531895!--       In case of non-flat topography define 2d-arrays containing the height
     
    21682310                CALL netcdf_handle_error( 'netcdf_define_header', 86 )
    21692311             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 )
    21702366
    21712367          ENDIF
     
    26732869                                 (/ id_dim_x_xy(av) /),      &
    26742870                                 'E_UTM', NF90_DOUBLE, id_var_eutm_xy(av,0),  &
    2675                                  'm', '', 000, 000, 000 )
     2871                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    26762872             CALL netcdf_create_var( id_set_xy(av), &
    26772873                                 (/ id_dim_y_xy(av) /),      &
    26782874                                 'N_UTM', NF90_DOUBLE, id_var_nutm_xy(av,0),  &
    2679                                  'm', '', 000, 000, 000 )
     2875                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    26802876             CALL netcdf_create_var( id_set_xy(av), &
    26812877                                 (/ id_dim_xu_xy(av) /),     &
    26822878                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_xy(av,1), &
    2683                                  'm', '', 000, 000, 000 )
     2879                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    26842880             CALL netcdf_create_var( id_set_xy(av), &
    26852881                                 (/ id_dim_y_xy(av) /),     &
    26862882                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_xy(av,1), &
    2687                                  'm', '', 000, 000, 000 )
     2883                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    26882884             CALL netcdf_create_var( id_set_xy(av), &
    26892885                                 (/ id_dim_x_xy(av) /),     &
    26902886                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_xy(av,2), &
    2691                                  'm', '', 000, 000, 000 )
     2887                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    26922888             CALL netcdf_create_var( id_set_xy(av), &
    26932889                                 (/ id_dim_yv_xy(av) /),     &
    26942890                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_xy(av,2), &
    2695                                  'm', '', 000, 000, 000 )
     2891                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    26962892          ELSE
    26972893             CALL netcdf_create_var( id_set_xy(av), &
    26982894                                 (/ id_dim_x_xy(av), id_dim_y_xy(av) /),      &
    26992895                                 'E_UTM', NF90_DOUBLE, id_var_eutm_xy(av,0),  &
    2700                                  'm', '', 000, 000, 000 )
     2896                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    27012897             CALL netcdf_create_var( id_set_xy(av), &
    27022898                                 (/ id_dim_x_xy(av), id_dim_y_xy(av) /),      &
    27032899                                 'N_UTM', NF90_DOUBLE, id_var_nutm_xy(av,0),  &
    2704                                  'm', '', 000, 000, 000 )
     2900                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    27052901             CALL netcdf_create_var( id_set_xy(av), &
    27062902                                 (/ id_dim_xu_xy(av), id_dim_y_xy(av) /),     &
    27072903                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_xy(av,1), &
    2708                                  'm', '', 000, 000, 000 )
     2904                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    27092905             CALL netcdf_create_var( id_set_xy(av), &
    27102906                                 (/ id_dim_xu_xy(av), id_dim_y_xy(av) /),     &
    27112907                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_xy(av,1), &
    2712                                  'm', '', 000, 000, 000 )
     2908                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    27132909             CALL netcdf_create_var( id_set_xy(av), &
    27142910                                 (/ id_dim_x_xy(av), id_dim_yv_xy(av) /),     &
    27152911                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_xy(av,2), &
    2716                                  'm', '', 000, 000, 000 )
     2912                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    27172913             CALL netcdf_create_var( id_set_xy(av), &
    27182914                                 (/ id_dim_x_xy(av), id_dim_yv_xy(av) /),     &
    27192915                                 '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 )
    27222947!
    27232948!--       In case of non-flat topography define 2d-arrays containing the height
     
    31993424          ENDIF
    32003425!
     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!
    32013480!--       In case of non-flat topography write height information. Only for
    32023481!--       parallel netcdf output.
     
    35953874                                 (/ id_dim_x_xz(av) /),      &
    35963875                                 'E_UTM', NF90_DOUBLE, id_var_eutm_xz(av,0),  &
    3597                                  'm', '', 000, 000, 000 )
     3876                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    35983877             CALL netcdf_create_var( id_set_xz(av), &
    35993878                                 (/ id_dim_y_xz(av) /),      &
    36003879                                 'N_UTM', NF90_DOUBLE, id_var_nutm_xz(av,0),  &
    3601                                  'm', '', 000, 000, 000 )
     3880                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    36023881             CALL netcdf_create_var( id_set_xz(av), &
    36033882                                 (/ id_dim_xu_xz(av) /),     &
    36043883                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_xz(av,1), &
    3605                                  'm', '', 000, 000, 000 )
     3884                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    36063885             CALL netcdf_create_var( id_set_xz(av), &
    36073886                                 (/ id_dim_y_xz(av) /),     &
    36083887                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_xz(av,1), &
    3609                                  'm', '', 000, 000, 000 )
     3888                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    36103889             CALL netcdf_create_var( id_set_xz(av), &
    36113890                                 (/ id_dim_x_xz(av) /),     &
    36123891                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_xz(av,2), &
    3613                                  'm', '', 000, 000, 000 )
     3892                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    36143893             CALL netcdf_create_var( id_set_xz(av), &
    36153894                                 (/ id_dim_yv_xz(av) /),     &
    36163895                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_xz(av,2), &
    3617                                  'm', '', 000, 000, 000 )
     3896                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    36183897          ELSE
    36193898             CALL netcdf_create_var( id_set_xz(av), &
    36203899                                 (/ id_dim_x_xz(av), id_dim_y_xz(av) /),      &
    36213900                                 'E_UTM', NF90_DOUBLE, id_var_eutm_xz(av,0),  &
    3622                                  'm', '', 000, 000, 000 )
     3901                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    36233902             CALL netcdf_create_var( id_set_xz(av), &
    36243903                                 (/ id_dim_x_xz(av), id_dim_y_xz(av) /),      &
    36253904                                 'N_UTM', NF90_DOUBLE, id_var_nutm_xz(av,0),  &
    3626                                  'm', '', 000, 000, 000 )
     3905                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    36273906             CALL netcdf_create_var( id_set_xz(av), &
    36283907                                 (/ id_dim_xu_xz(av), id_dim_y_xz(av) /),     &
    36293908                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_xz(av,1), &
    3630                                  'm', '', 000, 000, 000 )
     3909                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    36313910             CALL netcdf_create_var( id_set_xz(av), &
    36323911                                 (/ id_dim_xu_xz(av), id_dim_y_xz(av) /),     &
    36333912                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_xz(av,1), &
    3634                                  'm', '', 000, 000, 000 )
     3913                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    36353914             CALL netcdf_create_var( id_set_xz(av), &
    36363915                                 (/ id_dim_x_xz(av), id_dim_yv_xz(av) /),     &
    36373916                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_xz(av,2), &
    3638                                  'm', '', 000, 000, 000 )
     3917                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    36393918             CALL netcdf_create_var( id_set_xz(av), &
    36403919                                 (/ id_dim_x_xz(av), id_dim_yv_xz(av) /),     &
    36413920                                 '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 )
    36443952
    36453953          IF ( land_surface )  THEN
     
    40804388                DEALLOCATE( netcdf_data_2d )
    40814389             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 )
    40824445
    40834446          ENDIF
     
    44444807                                 (/ id_dim_x_yz(av) /),      &
    44454808                                 'E_UTM', NF90_DOUBLE, id_var_eutm_yz(av,0),  &
    4446                                  'm', '', 000, 000, 000 )
     4809                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    44474810             CALL netcdf_create_var( id_set_yz(av), &
    44484811                                 (/ id_dim_y_yz(av) /),      &
    44494812                                 'N_UTM', NF90_DOUBLE, id_var_nutm_yz(av,0),  &
    4450                                  'm', '', 000, 000, 000 )
     4813                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    44514814             CALL netcdf_create_var( id_set_yz(av), &
    44524815                                 (/ id_dim_xu_yz(av) /),     &
    44534816                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_yz(av,1), &
    4454                                  'm', '', 000, 000, 000 )
     4817                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    44554818             CALL netcdf_create_var( id_set_yz(av), &
    44564819                                 (/ id_dim_y_yz(av) /),     &
    44574820                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_yz(av,1), &
    4458                                  'm', '', 000, 000, 000 )
     4821                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    44594822             CALL netcdf_create_var( id_set_yz(av), &
    44604823                                 (/ id_dim_x_yz(av) /),     &
    44614824                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_yz(av,2), &
    4462                                  'm', '', 000, 000, 000 )
     4825                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    44634826             CALL netcdf_create_var( id_set_yz(av), &
    44644827                                 (/ id_dim_yv_yz(av) /),     &
    44654828                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_yz(av,2), &
    4466                                  'm', '', 000, 000, 000 )
     4829                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    44674830          ELSE
    44684831             CALL netcdf_create_var( id_set_yz(av), &
    44694832                                 (/ id_dim_x_yz(av), id_dim_y_yz(av) /),      &
    44704833                                 'E_UTM', NF90_DOUBLE, id_var_eutm_yz(av,0),  &
    4471                                  'm', '', 000, 000, 000 )
     4834                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    44724835             CALL netcdf_create_var( id_set_yz(av), &
    44734836                                 (/ id_dim_x_yz(av), id_dim_y_yz(av) /),      &
    44744837                                 'N_UTM', NF90_DOUBLE, id_var_nutm_yz(av,0),  &
    4475                                  'm', '', 000, 000, 000 )
     4838                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    44764839             CALL netcdf_create_var( id_set_yz(av), &
    44774840                                 (/ id_dim_xu_yz(av), id_dim_y_yz(av) /),     &
    44784841                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_yz(av,1), &
    4479                                  'm', '', 000, 000, 000 )
     4842                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    44804843             CALL netcdf_create_var( id_set_yz(av), &
    44814844                                 (/ id_dim_xu_yz(av), id_dim_y_yz(av) /),     &
    44824845                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_yz(av,1), &
    4483                                  'm', '', 000, 000, 000 )
     4846                                 'm', 'projection_y_coordinate', 000, 000, 000 )
    44844847             CALL netcdf_create_var( id_set_yz(av), &
    44854848                                 (/ id_dim_x_yz(av), id_dim_yv_yz(av) /),     &
    44864849                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_yz(av,2), &
    4487                                  'm', '', 000, 000, 000 )
     4850                                 'm', 'projection_x_coordinate', 000, 000, 000 )
    44884851             CALL netcdf_create_var( id_set_yz(av), &
    44894852                                 (/ id_dim_x_yz(av), id_dim_yv_yz(av) /),     &
    44904853                                 '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 )
    44934885
    44944886          IF ( land_surface )  THEN
     
    49165308                DEALLOCATE( netcdf_data_2d )
    49175309             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 )
    49185365
    49195366          ENDIF
     
    68967343 END SUBROUTINE netcdf_create_var
    68977344
     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
     7415END 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
    68987558 END MODULE netcdf_interface
Note: See TracChangeset for help on using the changeset viewer.