# **************************** LICENSE START ***********************************
#
# Copyright 2012 ECMWF. This software is distributed under the terms
# of the Apache License version 2.0. In applying this license, ECMWF does not
# waive the privileges and immunities granted to it by virtue of its status as
# an Intergovernmental Organization or submit itself to any jurisdiction.
#
# ***************************** LICENSE END ************************************

extern gradientb(f:fieldset) "fortran90" inline
!
!  "GRADIENTB" COMPUTES THE GRADIENT OF A FIELD
!
!  THIS PROGRAM IS A MODIFIED VERSION OF THE FORMER
!  "GRADIENT" TO TAKE INTO ACCOUNT THE UNITS.
! 
!  THE UNITS ARE IN THE INTERNATIONAL SYSTEM
!
!  GRIBEX version:    October,  1996
!  GRIB_API version:  March,    2010
!  MFI version:       November, 2010

  PROGRAM GRADIENTB
  USE grib_api
  IMPLICIT NONE

  INTEGER  fieldset_in, fieldset_out, icnt
  INTEGER  grib_id, isize, istatus, i
  INTEGER  byte_size
  REAL*8,  ALLOCATABLE :: grib_in(:)
  REAL*8,  ALLOCATABLE :: grib_out_u(:)
  REAL*8,  ALLOCATABLE :: grib_out_v(:)

                                    !-- GET FIRST ARGUMENT AS A FIELDSET.
                                    !-- icnt IS THE NUMBER OF FIELDS
  CALL mfi_get_fieldset( fieldset_in, icnt )

                                    !-- CREATE A NEW OUTPUT FIELDSET
  CALL mfi_new_fieldset( fieldset_out )

                                    !-- LOOP ON FIELDS
  DO i=1, icnt
                                    !-- GET NEXT FIELD FROM INPUT FIELDSET
     CALL mfi_load_one_grib( FIELDSET_IN, grib_id )

                                    !-- ALLOCATE ARRAYS, GET FIELD VALUES
     CALL grib_get_size( grib_id, 'values', isize )
     ALLOCATE( grib_in(isize), grib_out_u(isize), grib_out_v(isize) )

     CALL grib_get_real8_array( grib_id, 'values', grib_in, istatus )

                                    !-- VALIDATE AND DERIVE OUTPUT
     CALL valid( grib_id )
     CALL grad( grib_in, grib_out_u, grib_out_v )

                                    !-- SET OUTPUT AS U COMPONENT OF WIND
     CALL grib_set_real8_array( grib_id, 'values', grib_out_u, istatus )
     CALL grib_set_int( grib_id, 'paramId', 131 )

                                    !-- ADD IT TO THE OUTPUT FIELDSET
     CALL mfi_save_grib( fieldset_out, grib_id )

                                    !-- SET OUTPUT AS V COMPONENT OF WIND
     CALL grib_set_real8_array( grib_id, 'values', grib_out_v, istatus )
     CALL grib_set_int( grib_id, 'paramId', 132 )

                                    !-- ADD IT TO THE OUTPUT FIELDSET
     CALL mfi_save_grib( fieldset_out, grib_id )

                                    !-- RELEASE MEMORY
     CALL grib_release( grib_id )
     DEALLOCATE( grib_in, grib_out_u, grib_out_v )

  END DO
                                    !-- RETURN THE RESULT
  CALL mfi_return_fieldset( fieldset_out )

  STOP
END


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!--
!--  USER ROUTINE TO CHECK VALIDITY OF INPUT FIELD
!--  VALID FOR A GLOBAL FIELD, LAT/LONG, 1.5 DEG GRID
!--

SUBROUTINE valid( grib_id )
  USE grib_api
  INTEGER grib_id
  INTEGER ivalue
  REAL*8  rvalue

  CALL grib_get_int(grib_id, 'dataRepresentationType', ivalue)
  IF( ivalue .NE. 0 ) CALL mfi_fail("GRID not lat/lon")

  CALL grib_get_real8(grib_id, 'iDirectionIncrementInDegrees', rvalue)
  IF( rvalue .NE. 1.5 ) CALL mfi_fail("GRID not 1.5/1.5")

  CALL grib_get_real8(grib_id, 'jDirectionIncrementInDegrees', rvalue)
  IF( rvalue .NE. 1.5 ) CALL mfi_fail("GRID not 1.5/1.5")

  CALL grib_get_int( grib_id, 'Ni', ivalue )
  IF( ivalue .NE. 240 ) CALL mfi_fail("GRID not global")

  CALL grib_get_int( grib_id, 'Nj', ivalue )
  IF( ivalue .NE. 121 ) CALL mfi_fail("GRID not global")

  RETURN
END


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!--
!--  DERIVE GRADIENT OF INPUT FIELD F (BENITO ELVIRA, IM)
!--  FA = HORIZONTAL GRADIENT, FB = VERTICAL GRADIENT
!--

SUBROUTINE GRAD (F, FA, FB)

  !-- DIMENSIONS CORRESPONDING TO 1.5 x 1.5 GRID
  DIMENSION F(240,121), FA(240,121), FB(240,121)

  PI = ACOS(-1.0)
  RT = 6371000.0
  CB = (RT*1.5*PI)/180.0
                                    !-- COMPUTE HORIZONTAL GRADIENT
  DO I = 1, 121

     C = COS( (90.0-I*1.5 + 1.5)*PI/180.0 )
     FA(1,i) = (F(2,i)-F(240,i)) / (2.0*C*CB)
     FA(240,i) = (F(1,i)-F(239,i)) / (2.0*C*CB)

     DO J = 2, 239
        FA(j,i) = (F(j+1,i)-F(j-1,i)) / (2.0*C*CB)
     END DO

  END DO
                                    !-- COMPUTE VERTICAL GRADIENT
  DO I = 1, 240

     FB(i,1) = 0
     FB(i,121) = 0
     DO J = 2, 120
        FB(i,j) = (F(i,j+1)-F(i,j-1)) / (-2.0*CB)
     END DO

  END DO

  RETURN
END

end inline


