ol_write_proj.F90

Go to the documentation of this file.

SUBROUTINE ol_write_proj(HSELECT, KFILE_ID, UG)
!
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
!
USE modn_io_offline, ONLY : lwrite_coord
!
USE modd_csts, ONLY : xpi
!
USE modi_def_var_netcdf
USE modi_abor1_sfx
!
USE mode_gridtype_ign
USE mode_gridtype_conf_proj
USE mode_gridtype_lonlat_rot
USE mode_gridtype_lonlat_reg
USE mode_gridtype_lonlatval
USE mode_gridtype_gauss
USE mode_gridtype_cartesian
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE netcdf
!
IMPLICIT NONE
! 
 CHARACTER(LEN=*), DIMENSION(:), INTENT(IN) :: HSELECT
INTEGER,               INTENT(IN) :: KFILE_ID
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
!
 CHARACTER(LEN=100), DIMENSION(1) :: YATT_TITLE3
 CHARACTER(LEN=100), DIMENSION(11) :: YATT_TITLE2
 CHARACTER(LEN=100), DIMENSION(3) :: YATT_TITLE, YATT
 CHARACTER(LEN=100) :: YCOMMENT
REAL, DIMENSION(11) :: ZATT
REAL, DIMENSION(1,2) :: ZATT2D
!
REAL, DIMENSION(:), ALLOCATABLE :: ZLATP, ZLATP2
REAL, DIMENSION(:), ALLOCATABLE :: ZDX, ZDY, ZX, ZY
!
REAL, DIMENSION(2) :: ZLAT2
REAL, DIMENSION(1) :: ZLAT1, ZLON1, ZX1, ZY1, ZX2, ZY2
REAL :: ZLAT0, ZLON0, ZX0, ZY0, ZRPK, ZBETA, ZLATOR, ZLONOR, ZSIN, &
        ZLONCEN, ZLATCEN, ZCODIL, ZRESX, ZRESY, ZLONMIN, ZLONMAX,  &
        ZLATMIN, ZLATMAX, ZDX0, ZDY0
!
INTEGER, DIMENSION(:), ALLOCATABLE :: INLOPA
INTEGER, DIMENSION(1) :: IDIMID1
INTEGER, DIMENSION(0) :: IDIMID
INTEGER :: ILAMBERT, IL, IRESP, IVAR_ID, IIMAX, IJMAX, JJ, INLATI
INTEGER :: IIL, IIH, IJL, IJH, JRET, ICPT, ISIZE
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('OL_WRITE_PROJ',0,zhook_handle)
!
yatt_title(1) = "grid_mapping_name"
!
SELECT CASE(ug%G%CGRID)

  CASE ("CONF PROJ ")

    CALL get_gridtype_conf_proj(ug%XGRID_FULL_PAR,kl=il,kimax=iimax,kjmax=ijmax)
    ALLOCATE(zx(il),zy(il))
    ALLOCATE(zdx(il),zdy(il))
    CALL get_gridtype_conf_proj(ug%XGRID_FULL_PAR,zlat0,zlon0,zrpk,zbeta,zlator,zlonor,&
                                px=zx,py=zy,pdx=zdx,pdy=zdy)

    zresx = zdx(1)
    zresy = zdy(1)
    DEALLOCATE(zdx,zdy)

    IF ( abs(zrpk)==1 ) THEN
      yatt(1) = "polar_stereographic"
    ELSEIF (zrpk==0) THEN
      yatt(1) = "mercator"
    ELSE
      yatt(1) = "lambert_conformal_conic"
    ENDIF

    isize = 0

    ! geoid
    isize = isize + 1
    yatt_title2(isize) = 'earth_radius'
    zatt(isize) = 6371229.

    isize = isize + 1
    IF ( yatt(1)=="polar_stereographic" ) THEN
      yatt_title2(isize) = "straight_vertical_longitude_from_pole"
    ELSE
      yatt_title2(isize) = "longitude_of_central_meridian"
    ENDIF
    zatt(isize) = zlon0 

    isize = isize + 1
    yatt_title2(isize) = "latitude_of_projection_origin"
    IF (yatt(1)=="mercator") THEN
      zatt(isize) = 0.
    ELSE
      zatt(isize) = zlat0
    ENDIF

    zx1(1) = zx(1)
    zy1(1) = zy(1)    
    DEALLOCATE(zx,zy)    
    zlat1(1) = zlat0
    zlon1(1) = zlon0
    CALL xy_conf_proj(zlat0,zlon0,zrpk,0.,zlator,zlonor,zx2,zy2,zlat1,zlon1)
    isize = isize + 1
    yatt_title2(isize) = "false_easting" 
    zatt(isize) = nint((zx2(1)-zx1(1))*100.)/100.
    isize = isize + 1
    yatt_title2(isize) = "false_northing"
    zatt(isize) = nint((zy2(1)-zy1(1))*100.)/100.

    isize = isize + 1
    yatt_title2(isize) = "rotation"
    zatt(isize) = zbeta
    
    IF (zresx/=0.) THEN
      isize = isize + 1
      yatt_title2(isize) = "x_resolution"
      zatt(isize) = zresx
      isize = isize + 1
      yatt_title2(isize) = "y_resolution"
      zatt(isize) = zresy
    ENDIF
    
    zsin = sin(zlat0*(xpi/180.))
    IF (abs(zrpk)==1.OR.zrpk==0..OR.nint(zsin*100.)==nint(zrpk*100.)) THEN
    
      isize = isize + 1
      yatt_title2(isize) = 'standard_parallel'
      zatt(isize) = zlat0

      CALL def_var_netcdf(hselect,kfile_id,"Projection_Type","",idimid,&
            yatt_title(1:1),yatt(1:1),ivar_id,nf90_int,hatt_title2=yatt_title2(1:isize),&
            patt_float=zatt(1:isize))
      
    ELSE
       
      yatt_title3(1) = "standard_parallel"
      zatt2d(1,1) = asin(zrpk)*(180./xpi)
      zatt2d(1,2) = 2*zlat0-asin(zrpk)*(180./xpi)

      CALL def_var_netcdf(hselect,kfile_id,"Projection_Type","",idimid,&
          yatt_title(1:1),yatt(1:1),ivar_id,nf90_int,hatt_title2=yatt_title2(1:isize),&
          patt_float=zatt(1:isize),hatt_title3=yatt_title3(1:1),patt_float2d=zatt2d(1:1,:))
        
    ENDIF


  CASE ("GAUSS     ")

    yatt(1) = "gauss_grid"

    CALL get_gridtype_gauss(ug%XGRID_FULL_PAR,plapo=zlat0,plopo=zlon0,kl=il,knlati=inlati)
    ALLOCATE(inlopa(inlati))
    ALLOCATE(zlatp(il))
    ALLOCATE(zlatp2(inlati))
    CALL get_gridtype_gauss(ug%XGRID_FULL_PAR,knlopa=inlopa,pcodil=zcodil,plat_xy=zlatp)
    zlatp2(1) = zlatp(1)
    icpt = 1
    DO jj = 2,il
      IF ( zlatp(jj)/=zlatp(jj-1) ) THEN
        icpt = icpt + 1
        zlatp2(icpt) = zlatp(jj)
      ENDIF
    ENDDO
    DEALLOCATE(zlatp)

    yatt_title(2) = "latitudes"
    yatt_title(3) = "lon_number_by_lat"
    yatt(2) = "var: gauss_latitudes"
    yatt(3) = "var: lon_number_by_lat"
    
    ! geoid
    yatt_title2(1) = 'earth_radius'
    zatt(1) = 6371229.

    yatt_title2(2) = "dilatation_coef"
    zatt(2) = zcodil

    yatt_title2(3) = "pole_lon"
    yatt_title2(4) = "pole_lat"
    zatt(3) = nint(zlon0*100.)/100.
    zatt(4) = nint(zlat0*100.)/100.
    
    CALL def_var_netcdf(hselect,kfile_id,"Projection_Type","",idimid,&
              yatt_title(1:3),yatt(1:3),ivar_id,nf90_int,hatt_title2=yatt_title2(1:4),&
              patt_float=zatt(1:4))

    yatt_title(:) = ""
    yatt(:) = ""
    jret = nf90_inq_dimid(kfile_id,"latitude",idimid1(1))
    CALL def_var_netcdf(hselect,kfile_id,"gauss_latitudes","",idimid1,&
              yatt_title(1:1),yatt(1:1),ivar_id,nf90_double)
    CALL def_var_netcdf(hselect,kfile_id,"lon_number_by_lat","",idimid1,&
              yatt_title(1:1),yatt(1:1),ivar_id,nf90_int)
    jret = nf90_enddef(kfile_id)
    jret = nf90_inq_varid(kfile_id, "gauss_latitudes", ivar_id)
    jret = nf90_put_var(kfile_id, ivar_id, zlatp2)
    jret = nf90_inq_varid(kfile_id, "lon_number_by_lat", ivar_id)
    jret = nf90_put_var(kfile_id, ivar_id, inlopa)
    jret = nf90_redef(kfile_id)

    DEALLOCATE(zlatp2,inlopa)


  CASE ("CARTESIAN ")

    yatt(1) = "cartesian"

    CALL get_gridtype_cartesian(ug%XGRID_FULL_PAR,kl=il)
    ALLOCATE(zdx(il),zdy(il))
    CALL get_gridtype_cartesian(ug%XGRID_FULL_PAR,pdx=zdx,pdy=zdy)    
    zresx = zdx(1)
    zresy = zdy(1)
    DEALLOCATE(zdx,zdy)

    isize = 0
    IF (zresx/=0.) THEN
      isize = isize + 1
      yatt_title2(isize) = "x_resolution"
      zatt(isize) = zresx
      isize = isize + 1
      yatt_title2(isize) = "y_resolution"
      zatt(isize) = zresy
    ENDIF
!
    CALL def_var_netcdf(hselect,kfile_id,"Projection_Type","",idimid,&
                 yatt_title(1:1),yatt(1:1),ivar_id,nf90_int,hatt_title2=yatt_title2(1:isize),&
                 patt_float=zatt(1:isize))


  CASE ("LONLAT REG")

    yatt(1) = "latitude_longitude"

    CALL get_gridtype_lonlat_reg(ug%XGRID_FULL_PAR,kl=il,plonmin=zlonmin,&
                                 plonmax=zlonmax,platmin=zlatmin,platmax=zlatmax)
    ALLOCATE(zx(il),zy(il))
    CALL get_gridtype_lonlat_reg(ug%XGRID_FULL_PAR,plon=zx,plat=zy)

    IF (il>1) THEN
      zresx = abs(zx(1)-zx(2))
      iimax = nint((zlonmax-zlonmin)/zresx)
      IF (il>iimax) THEN
        zresy = abs(zy(1)-zy(iimax+1))
      ELSE
        zresy = zlatmax - zlatmin
      ENDIF
    ELSE
      zresx = zlonmax - zlonmin
      zresy = zlatmax - zlatmin
    ENDIF
    DEALLOCATE(zx,zy)

    ! geoid
    yatt_title2(1) = 'earth_radius'
    zatt(1) = 6371229.
    
    isize = 1
    IF (zresx/=0.) THEN
      isize = isize + 1
      yatt_title2(isize) = "x_resolution"
      zatt(isize) = zresx
      isize = isize + 1
      yatt_title2(isize) = "y_resolution"
      zatt(isize) = zresy
    ENDIF
    
    CALL def_var_netcdf(hselect,kfile_id,"Projection_Type","",idimid,&
                 yatt_title(1:1),yatt(1:1),ivar_id,nf90_int,hatt_title2=yatt_title2(1:isize),&     
                 patt_float=zatt(1:isize))


  CASE ("LONLAT ROT")

    yatt(1) = "rotated_latitude_longitude"

    CALL get_gridtype_lonlat_rot(ug%XGRID_FULL_PAR,ppolon=zlon0,ppolat=zlat0,&
                                 kl=il,klon=iimax,klat=ijmax,pdlon=zdx0,pdlat=zdy0)

    ! geoid
    yatt_title2(1) = 'earth_radius'
    zatt(1) = 6371229.
                         
    yatt_title2(2) = "grid_north_pole_longitude"    
    yatt_title2(3) = "grid_north_pole_latitude"
    zatt(2) = zlon0
    zatt(3) = zlat0

    zresx = zdx0
    zresy = zdy0

    isize = 3
    IF (zresx/=0.) THEN
      isize = isize + 1
      yatt_title2(isize) = "x_resolution"
      zatt(isize) = zresx
      isize = isize + 1
      yatt_title2(isize) = "y_resolution"
      zatt(isize) = zresy
    ENDIF
 
    CALL def_var_netcdf(hselect,kfile_id,"Projection_Type","",idimid,&
                  yatt_title(1:1),yatt(1:1),ivar_id,nf90_int,&
                  hatt_title2=yatt_title2(1:isize),patt_float=zatt(1:isize))
        

  CASE ("LONLATVAL ")

    yatt(1) = "latitude_longitude"

    ! geoid
    yatt_title2(1) = 'earth_radius'
    zatt(1) = 6371229.
    
    CALL def_var_netcdf(hselect,kfile_id,"Projection_Type","",idimid,&
                        yatt_title(1:1),yatt(1:1),ivar_id,nf90_int,&
                  hatt_title2=yatt_title2(1:1),patt_float=zatt(1:1))


  CASE ("IGN       ")

    CALL get_gridtype_ign(ug%XGRID_FULL_PAR,klambert=ilambert,kl=il,&
                         kdimx=iimax,kdimy=ijmax)

    IF ( lwrite_coord .AND. iimax*ijmax/=il ) THEN

      yatt(1) = "latitude_longitude"

      CALL def_var_netcdf(hselect,kfile_id,"Projection_Type","",idimid,&
                 yatt_title(1:1),yatt(1:1),ivar_id,nf90_int)

    ELSE

      ALLOCATE(zdx(il))
      ALLOCATE(zdy(il))
      ALLOCATE(zx(iimax),zy(ijmax))
      CALL get_gridtype_ign(ug%XGRID_FULL_PAR,pxall=zx,pyall=zy,pdx=zdx,pdy=zdy)

      SELECT CASE(ilambert)
        CASE(1)
          zlat0 = 49.5
          zlat2(1) = 48.592447
          zlat2(2) = 50.39088033333
          zlon0 = 2.33722917
          zx0 = 600000.
          zy0 = 5657616.674
            
        CASE(2)
          zlat0 = 46.8
          zlat2(1) = 45.88935133333
          zlat2(2) = 47.68760866666     
          zlon0 = 2.33722917
          zx0 = 600000.
          zy0 = 6199695.768

        CASE(3)
          zlat0 = 44.1
          zlat2(1) = 43.19290816666
          zlat2(2) = 44.96513966666    
          zlon0 = 2.33722917
          zx0 = 600000.
          zy0 = 6791905.085
          
        CASE(4)
          zlat0 = 42.165
          zlat2(1) = 41.55623266666
          zlat2(2) = 42.76726466666      
          zlon0 = 2.33722917
          zx0 = 234.358
          zy0 = 7239161.542
  
        CASE(5)
          zlat0 = 46.8
          zlat2(1) = 45.88935133333
          zlat2(2) = 47.68760866666     
          zlon0 = 2.33722917
          zx0 = 600000.
          zy0 = 8199695.768
  
        CASE(6)
          zlat0 = 46.5
          zlat2(1) = 44.
          zlat2(2) = 49.         
          zlon0 = 3.
          zx0 = 700000.
          zy0 = 12655612.050
            
      END SELECT

      yatt_title(1) = "grid_mapping_name"
      yatt(1) = "lambert_conformal_conic"
      yatt_title(2) = "ellipsoid"
      IF (ilambert==6) THEN
        yatt(2) = "GRS80"
      ELSE
        yatt(2) = "clrk80"
      ENDIF

      isize = 0

      isize = isize + 1
      yatt_title2(isize) = "longitude_of_central_meridian"
      zatt(isize) = zlon0 

      isize = isize + 1
      yatt_title2(isize) = "latitude_of_projection_origin"
      zatt(isize) = zlat0

      zlat1(1) = zlat0
      zlon1(1) = zlon0
      CALL xy_ign(ilambert,zx2,zy2,zlat1,zlon1)
      zx1(1) = minval(zx)
      zy1(1) = minval(zy)
      isize = isize + 1
      yatt_title2(isize) = "false_easting" 
      zatt(isize) = zx2(1)-zx1(1)
      isize = isize + 1
      yatt_title2(isize) = "false_northing"
      zatt(isize) = zy2(1)-zy1(1)
        
      yatt_title3(1) = "standard_parallel"
      zatt2d(1,:) = zlat2(:)
         
      IF (all(zdx(2:)==zdx(1))) THEN
        isize = isize + 1
        yatt_title2(isize) = "x_resolution"
        zatt(isize) = zdx(1)
      ENDIF
      IF (all(zdy(2:)==zdy(1))) THEN
        isize = isize + 1
        yatt_title2(isize) = "y_resolution"
        zatt(isize) = zdy(1)
      ENDIF

      CALL def_var_netcdf(hselect,kfile_id,"Projection_Type","",idimid,&
            yatt_title(1:2),yatt(1:2),ivar_id,nf90_int,hatt_title2=yatt_title2(1:isize),&
            patt_float=zatt(1:isize),hatt_title3=yatt_title3(1:1),patt_float2d=zatt2d(1:1,:))
  
      DEALLOCATE(zx,zy,zdx,zdy)

   ENDIF

  CASE DEFAULT

    CALL abor1_sfx("OL_WRITE_PROJ: PROJECTION "//ug%G%CGRID//" NOT DEFINED")

END SELECT 
!
IF (lhook) CALL dr_hook('OL_WRITE_PROJ',1,zhook_handle)
!
END SUBROUTINE ol_write_proj