MODULE CPLNG_DATA_MOD

    USE PARKIND1, ONLY: JPIM, JPRB

    USE CPLNG_LOG_MOD

    IMPLICIT NONE

    PRIVATE

    PUBLIC CPLNG_ACTIVE
    PUBLIC CPLNG_IS_ACTIVE

    PUBLIC CPLNG_FLD_TYPE
    PUBLIC CPLNG_FLD

    PUBLIC CPLNG_ADD_FLD
    PUBLIC CPLNG_ADD_FLD_COMPLETED
    PUBLIC CPLNG_IDX

    PUBLIC CPLNG_FLD_TYPE_GRIDPOINT
    PUBLIC CPLNG_FLD_TYPE_SPECTRAL

    LOGICAL :: CPLNG_ACTIVE = .FALSE.

    TYPE CPLNG_FLD_TYPE
        INTEGER(KIND=JPIM),ALLOCATABLE :: ID(:,:)
        CHARACTER(LEN=128)             :: NAME
        INTEGER(KIND=JPIM)             :: TYPE
        INTEGER(KIND=JPIM)             :: INOUT
        INTEGER(KIND=JPIM)             :: STAGE
        INTEGER(KIND=JPIM)             :: NUM_LVL
        INTEGER(KIND=JPIM)             :: NUM_CAT
        REAL(KIND=JPRB),   ALLOCATABLE :: D(:,:,:)
    END TYPE CPLNG_FLD_TYPE

    TYPE(CPLNG_FLD_TYPE),ALLOCATABLE,TARGET :: CPLNG_FLD(:)

    INTEGER(KIND=JPIM),PARAMETER :: CPLNG_FLD_TYPE_GRIDPOINT = 0
    INTEGER(KIND=JPIM),PARAMETER :: CPLNG_FLD_TYPE_SPECTRAL  = 1

    ! NUM_FIELDS is for internal bookkeeping
    INTEGER(KIND=JPIM) :: NUM_FIELDS = 0

CONTAINS

FUNCTION CPLNG_IS_ACTIVE

    LOGICAL :: CPLNG_IS_ACTIVE

    CPLNG_IS_ACTIVE = CPLNG_ACTIVE

END FUNCTION CPLNG_IS_ACTIVE

SUBROUTINE CPLNG_ADD_FLD(NAME,TYPE,INOUT,STAGE,LVL,CAT)

    ! PARAMETERS 
    !Number of coupling fields. Should be namelist variable?
    INTEGER,PARAMETER :: ALLOCATE_CHUNK = 20  !!Joakim edit 

    ! ARGUMENTS
    CHARACTER(LEN=*),  INTENT(IN)          :: NAME
    INTEGER(KIND=JPIM),INTENT(IN)          :: TYPE
    INTEGER(KIND=JPIM),INTENT(IN)          :: INOUT
    INTEGER(KIND=JPIM),INTENT(IN)          :: STAGE
    INTEGER(KIND=JPIM),INTENT(IN),OPTIONAL :: LVL
    INTEGER(KIND=JPIM),INTENT(IN),OPTIONAL :: CAT

    ! LOCALS
    INTEGER(KIND=JPIM)               :: NUM_LVL,NUM_CAT
    TYPE(CPLNG_FLD_TYPE),ALLOCATABLE :: CPLNG_FLD_TMP(:)

    ! Early return if CPLNG not active
    IF (.NOT.CPLNG_ACTIVE) THEN
        CALL CPLNG_WARNING('CPLNG_ADD_FLD called but CPLNG_ACTIVE is false.')
        RETURN
    ENDIF

    ! Set number of levels/categories to default or optional arguments
    NUM_LVL = 1
    NUM_CAT = 1

    IF (PRESENT(LVL)) NUM_LVL = LVL
    IF (PRESENT(CAT)) NUM_CAT = CAT

    ! Make initial allocation, if necessary
    IF (.NOT.ALLOCATED(CPLNG_FLD)) THEN
        ALLOCATE(CPLNG_FLD(ALLOCATE_CHUNK))
    ENDIF

    ! Increase (by reallocation) size of CPLNG_FLD, if necessary
    IF (NUM_FIELDS==SIZE(CPLNG_FLD)) THEN
        ALLOCATE(CPLNG_FLD_TMP(SIZE(CPLNG_FLD)+ALLOCATE_CHUNK))
        CPLNG_FLD_TMP(1:NUM_FIELDS) = CPLNG_FLD
        CALL MOVE_ALLOC(CPLNG_FLD_TMP,CPLNG_FLD) ! Fortran 2003
    ENDIF

    NUM_FIELDS=NUM_FIELDS+1

    CPLNG_FLD(NUM_FIELDS)%NAME    = NAME
    CPLNG_FLD(NUM_FIELDS)%TYPE    = TYPE
    CPLNG_FLD(NUM_FIELDS)%INOUT   = INOUT
    CPLNG_FLD(NUM_FIELDS)%STAGE   = STAGE
    CPLNG_FLD(NUM_FIELDS)%NUM_LVL = NUM_LVL
    CPLNG_FLD(NUM_FIELDS)%NUM_CAT = NUM_CAT

END SUBROUTINE CPLNG_ADD_FLD

FUNCTION CPLNG_IDX(FLD_NAME)
    ! Argument
    CHARACTER(LEN=*), INTENT(IN) :: FLD_NAME
    ! Result
    INTEGER(KIND=JPIM) :: CPLNG_IDX
    ! Locals
    INTEGER(KIND=JPIM) :: II

    ! Early return if CPLNG not active
    IF (.NOT.CPLNG_ACTIVE) THEN
        CALL CPLNG_WARNING('CPLNG_IDX called but CPLNG_ACTIVE is false.')
        RETURN
    ENDIF

    IF (.NOT.ALLOCATED(CPLNG_FLD)) THEN
        CALL ABOR1("CPLNG_IDX: CPLNG_FLD not allocated upon call")
    ENDIF

    ! Simple loop searching for the right field name
    DO II=1,SIZE(CPLNG_FLD)
        IF (TRIM(FLD_NAME) == TRIM(CPLNG_FLD(II)%NAME)) EXIT
    ENDDO

    IF (II>SIZE(CPLNG_FLD)) THEN
        CALL ABOR1("CPLNG_IDX: Field name not found: "//TRIM(FLD_NAME))
    ENDIF

    CPLNG_IDX = II
END FUNCTION CPLNG_IDX

SUBROUTINE CPLNG_ADD_FLD_COMPLETED
    
    ! NULOUT: id of ifs.log file
    ! NULNAM: id of namelist file
    USE YOMLUN, ONLY: NULOUT, NULNAM
    ! NGPTOT: size of arrays in this process (MPI task)
    ! NLOENG: number of longitudes on each latitude row
    USE YOMGEM, ONLY: NGPTOT, NGPTOTG, NLOENG 
    USE YOMLEG, ONLY: YRCSGLEG
    ! Constants (pi, radius of Earth)
    USE YOMCST, ONLY: RPI, RA
    USE YOMDIM, ONLY: NSMAX, NUMP, NSPEC2, NPROMA, NDGLG
    ! mapping from local to global indices
    USE YOMMP,  ONLY: NGLOBALINDEX 
    USE YOMMP0, ONLY: MYPROC
    USE YOMCT0, ONLY: NPROC
    USE YOMGC,  ONLY: YRGSGEOM_NB ! contains grid info
    USE SPECTRAL_FIELDS_DATA, ONLY: MYMS
    USE SURFACE_FIELDS_MIX, ONLY: SP_SB, YSP_SB, SP_RR, YSP_RR, SD_VF, YSD_VF
    
    USE MOD_OASIS

    ! Variables needed in OASIS calls
    INTEGER(KIND=JPIM)              :: IERROR, IERR, IOASIS3LVL
    INTEGER(KIND=JPIM), ALLOCATABLE :: IG_PARAL(:)
    INTEGER(KIND=JPIM)              :: IL_PART_ID, IL_PART_ID_GP, IL_PART_ID_SP
    INTEGER(KIND=JPIM)              :: IVAR_NODIMS(2)
    INTEGER(KIND=JPIM)              :: IVAR_ACTUAL_SHAPE(2)
    CHARACTER(LEN=4)                :: CGRDNAM = 'atmo'
    
    ! These two should be namelist variables or environment variables
    LOGICAL            :: LL_DEBUGCPL=.TRUE.
    LOGICAL            :: LL_SKIPGRIDWRITE=.FALSE.
    REAL(KIND=JPRB)    :: RLSMCRIT = 0.5_JPRB
    REAL(KIND=JPRB)    :: RNEMOSOUTH=-77.0105_JPRB
     
    LOGICAL            :: LSEGSTART(NGPTOT)
    INTEGER(KIND=JPIM) :: II,IP,IJ,JJ,JI,IJG,JKGLO,IST,IEND,I,J,IBLK
    INTEGER(KIND=JPIM) :: ILVL,ICAT
    INTEGER(KIND=JPIM) :: IOFFSET
    INTEGER(KIND=JPIM) :: IEXTENT
    INTEGER(KIND=JPIM) :: INUMSEGMENTS
    INTEGER(KIND=JPIM) :: IPOINT, ISTLAT, IENDLAT, ISTLON, IENDLON
    CHARACTER(LEN=128) :: CFLDNAME
    CHARACTER(LEN=3)   :: CERRSTR
    
    REAL(KIND=JPRB)                 :: ZDELTAX, ZDELTAYUP, ZDELTAYDW, ZGLOAREA, ZGLOAREA2
    REAL(KIND=JPRB),    ALLOCATABLE :: ZLON(:,:)   ! local longitudes
    REAL(KIND=JPRB),    ALLOCATABLE :: ZLAT(:,:)   ! local latitudes
    REAL(KIND=JPRB),    ALLOCATABLE :: ZLSM(:,:)   ! local land-sea mask (in fractions) 
    REAL(KIND=JPRB),    ALLOCATABLE :: ZAREA(:,:)  ! local grid-cell area
    REAL(KIND=JPRB),    ALLOCATABLE :: ZAREA2(:,:) 
    REAL(KIND=JPRB),    ALLOCATABLE :: ZCLAT(:,:,:)! latitudes of corner points for each cell  
    REAL(KIND=JPRB),    ALLOCATABLE :: ZCLON(:,:,:)! longitudes of corner points for each cell  
    REAL(KIND=JPRB),    ALLOCATABLE :: ZRGAUSLAT(:)! latitudes for each latitude row 
    INTEGER(KIND=JPIM), ALLOCATABLE :: ILSM(:,:)   ! local land-sea mask (binary, 0 or 1) 
    INTEGER(KIND=JPIM), ALLOCATABLE :: ILSM2(:,:)   ! local land-sea mask (binary, 0 or 1)
    
    INTEGER(KIND=JPIM) :: ISPM
    
#include "posnam.intfb.h"
#include "namcplng.nam.h" 
    
    CALL POSNAM(NULNAM,'NAMCPLNG') 
    READ       (NULNAM, NAMCPLNG)  
    
    IF (LL_DEBUGCPL) THEN
        CALL CPLNG_INFO('OpenIFS WILL NOW RUN THE OASIS COMMANDS TO DEFINE COUPLING FIELDS ')
    END IF
    ! Early return if CPLNG not active
    IF (.NOT.CPLNG_ACTIVE) THEN
        CALL CPLNG_WARNING('CPLNG_ADD_FLD_COMPLETED called but CPLNG_ACTIVE is false.')
        RETURN
    ENDIF

    ! -------------------------------------------------------------------------
    ! * (3) SET UP OASIS PARTITION
    ! -------------------------------------------------------------------------

    ! * (3A) Grid point fields
    IF( ANY(CPLNG_FLD(:)%TYPE == CPLNG_FLD_TYPE_GRIDPOINT) )THEN
        
        IF (LL_DEBUGCPL) THEN
            CALL CPLNG_INFO('OpenIFS: AT LEAST ONE COUPLING FIELD IS GRIDPOINT. ')
        END IF
        ! Set up logical array LSEGSTART that indicats where OASIS segments start
        ! Method: Use the NGLOBALINDEX array from YOMMP that contains a
        !         local-to-global mapping of the grid points. If two consecutive
        !         indices in this array differ by any other value than one, it's the
        !         start of a new segment.
        LSEGSTART(1)        = .TRUE.
        LSEGSTART(2:NGPTOT) = NGLOBALINDEX(2:NGPTOT)-NGLOBALINDEX(1:NGPTOT-1) /= 1

        INUMSEGMENTS = COUNT(LSEGSTART) ! Counts number of segments

        ! For shape and meaning of IG_PARAL, see OASIS documentation
        ALLOCATE(IG_PARAL(2 + 2*INUMSEGMENTS))

        IG_PARAL(1) = 3 ! Value 3 indicates an orange partition, which is
                        ! an ensemble of segments of the global domain
        IG_PARAL(2) = INUMSEGMENTS

        ! Compute segment offsets and extents and store them in IG_PARAL for later
        ! use in the OASIS_DEF_PARTITION call
        IOFFSET = NGLOBALINDEX(1)
        IEXTENT = 1
        IP      = 3 ! Pointer into IG_PARAL array
        DO II=2,NGPTOT
            IF (LSEGSTART(II)) THEN
                ! OASIS counts offsets starting with zero, hence minus one
                IG_PARAL(IP)   = IOFFSET-1
                IG_PARAL(IP+1) = IEXTENT
                IP = IP+2
                IOFFSET = NGLOBALINDEX(II)
                IEXTENT = 1
            ELSE
                IEXTENT = IEXTENT+1
            ENDIF
        ENDDO
        IG_PARAL(IP)   = IOFFSET-1
        IG_PARAL(IP+1) = IEXTENT

        ! Define partition for OASIS
        IF (LL_DEBUGCPL) THEN
            CALL CPLNG_INFO('OpenIFS: Define partition for gridpoint variable OASIS')
        END IF
        CALL OASIS_DEF_PARTITION(IL_PART_ID_GP,IG_PARAL,IERROR)
        IF (IERROR/=OASIS_OK) THEN
            WRITE (CERRSTR,'(I3)') IERROR
            CALL ABOR1("CPLNG_OCE_CONFIG: Error on OASIS_DEF_PARTITION (gridpoint): "//CERRSTR)
        ENDIF

        DEALLOCATE(IG_PARAL)

    ENDIF ! Field type == gridpoint

    ! * (3B) Spectral fields
    IF( ANY(CPLNG_FLD(:)%TYPE == CPLNG_FLD_TYPE_SPECTRAL) )THEN
       !
       ! Orange partition: each processor owns a number of (not following) wave numbers.
       !
       !  yomdim/NSMAX        : Spectral truncation T
       !  yomdim/NUMP         : Number of wavenumbers locally
       !  yomlap/MYMS(1:NUMP) : Actual wave numbers   (NOTE: not 0:NUMP as in comment!)
       
       IF (LL_DEBUGCPL) THEN
           CALL CPLNG_INFO('OpenIFS: At least one of coupling fields is spectral. Allocate space')
       END IF
       ! Allocate definition array
       ALLOCATE( IG_PARAL(2+2*NUMP) )

       ! Indicator for orange partition
       IG_PARAL(1) = 3

       ! Number of segments
       IG_PARAL(2) = NUMP

       ! Global offset and size
       DO II = 1, NUMP
          ! Global wave number
          ISPM = MYMS(II)
          ! Offset
          IG_PARAL(1+II*2) = ((NSMAX+1)+(NSMAX+2-ISPM)) * ISPM
          ! Size
          IG_PARAL(2+II*2) = (NSMAX+1-ISPM)*2  ! M:T, Re/Im
       END DO

       ! Define partition for OASIS
       IF (LL_DEBUGCPL) THEN
           CALL CPLNG_INFO('OpenIFS: Define partition in OASIS for spectral variable')
       END IF
       CALL OASIS_DEF_PARTITION(IL_PART_ID_SP,IG_PARAL,IERROR)
       IF (IERROR/=OASIS_OK) THEN
          WRITE (CERRSTR,'(I3)') IERROR
          CALL ABOR1("ECE_CPLNG_OCE_CONFIG: Error on OASIS_DEF_PARTITION (spectral): "//CERRSTR)
       ENDIF

       DEALLOCATE(IG_PARAL)

    ENDIF ! Field type == spectral

    ! -------------------------------------------------------------------------
    ! * (4) DEFINE COUPLING FIELDS FOR OASIS
    ! -------------------------------------------------------------------------
    IVAR_NODIMS = (/ 1, 1 /)

    WRITE (NULOUT,'(/,1X,A,/,1X,A,I4)')        &
    &     'CPLNG: Coupling field definition.', &
    &     'CPLNG: ... Number of fields is: ',SIZE(CPLNG_FLD)

    DO II=1,SIZE(CPLNG_FLD)

        WRITE (NULOUT,'(1X,A,I4,A,I2,A,I2,A,A20)') &
        &     'CPLNG: .... Field ',II,             &
        &     '(',CPLNG_FLD(II)%NUM_LVL,' lvl,',   &
        &     CPLNG_FLD(II)%NUM_CAT,' cat) Name: ',TRIM(CPLNG_FLD(II)%NAME)

        ! Allocate for ID's of all levels/categories of the field
        ALLOCATE(CPLNG_FLD(II)%ID(CPLNG_FLD(II)%NUM_LVL,CPLNG_FLD(II)%NUM_CAT))

        ! Set OASIS partition id and data size according to field type
        ! (gp/spectral). Allocate data member.
        SELECT CASE (CPLNG_FLD(II)%TYPE)

        CASE (CPLNG_FLD_TYPE_GRIDPOINT)
            IL_PART_ID        = IL_PART_ID_GP
            IVAR_ACTUAL_SHAPE = (/ 1, NGPTOT /)
            ALLOCATE(CPLNG_FLD(II)%D(NGPTOT,                &
            &                        CPLNG_FLD(II)%NUM_LVL, &
            &                        CPLNG_FLD(II)%NUM_CAT  ))

        CASE (CPLNG_FLD_TYPE_SPECTRAL )
            IL_PART_ID        = IL_PART_ID_SP
            IVAR_ACTUAL_SHAPE = (/ 1, NSPEC2 /)
            ALLOCATE(CPLNG_FLD(II)%D(NSPEC2,                &
            &                        CPLNG_FLD(II)%NUM_LVL, &
            &                        CPLNG_FLD(II)%NUM_CAT  ))

        CASE DEFAULT
            CALL ABOR1("CPLNG_OCE_CONFIG: Wrong field type for " // TRIM(CPLNG_FLD(II)%NAME))

        END SELECT

        ! Initialise to a strange value
        CPLNG_FLD(II)%D = -HUGE(0.0_JPRB)

        DO ICAT=1,CPLNG_FLD(II)%NUM_CAT
            DO ILVL=1,CPLNG_FLD(II)%NUM_LVL

                CFLDNAME = TRIM(CPLNG_FLD(II)%NAME)
                IF (CPLNG_FLD(II)%NUM_CAT>1) WRITE (CFLDNAME,'(A,".C",I3.3)') TRIM(CFLDNAME),ICAT
                IF (CPLNG_FLD(II)%NUM_LVL>1) WRITE (CFLDNAME,'(A,".L",I3.3)') TRIM(CFLDNAME),ILVL
                
                IF (LL_DEBUGCPL) THEN
                    CALL CPLNG_INFO('OpenIFS: Running OASIS_DEF_VAR for each coupling field '&
                                   & //TRIM(CFLDNAME))
                END IF
                CALL OASIS_DEF_VAR(CPLNG_FLD(II)%ID(ILVL,ICAT), &
                &                  CFLDNAME,                    &
                &                  IL_PART_ID,                  &
                &                  IVAR_NODIMS,                 &
                &                  CPLNG_FLD(II)%INOUT,         &
                &                  IVAR_ACTUAL_SHAPE,           &
                &                  OASIS_REAL,                  &
                &                  IERROR                       )
                WRITE (NULOUT,'(1X,A,I4,A,I3)') &
                &     'CPLNG: .... Field ',II,  &
                &     ' OASIS_DEF_VAR returns ',IERROR
                
                IF (IERROR/=OASIS_OK) THEN
                    WRITE (CERRSTR,'(I3)') IERROR
                    CALL ABOR1("CPLNG_OCE_CONFIG: Error in OASIS_DEF_VAR: " // CERRSTR // &
                    &          " (" // TRIM(CFLDNAME) // ")")
                ENDIF
            ENDDO
        ENDDO
    ENDDO
    
    ! -------------------------------------------------------------------------
    ! * (5) WRITE GRID INFORMATION, IF NEEDED
    ! Code taken from OpenIFS CY43R3, but modified for OpenIFS CY40R1
    ! -------------------------------------------------------------------------
    
    ! Check if we need to write grid information
    
    IERR=0
    IOASIS3LVL = 1 ! In CY43 we can set OASIS3LVL=0 to avoid using OASIS3
                   ! Here we let it always be = 1
    
    IF (.NOT. LL_SKIPGRIDWRITE) THEN
    
    !IF (MYPROC == 1) THEN
    
    IF (IOASIS3LVL == 1) THEN        
        CALL OASIS_START_GRIDS_WRITING(IERR) 
        ! In the OASIS3 code, grids writing ALWAYS returns IERR=1 so that
        ! grids are always written. 
        ! But we can only do it if we are on single-processor runs. 
        ! So when grids have been written, set LL_SKIPGRIDWRITE = .true.
        ! to skip grids writing. 
        !IF (LL_SKIPGRIDWRITE) THEN
        !   !IERR=0
        !END IF
        IF (LL_DEBUGCPL > 0) THEN
            WRITE(NULOUT,'(A,I4)') &
                 & 'Return code from oasis_start_grids_writing : ',IERR
        ENDIF
        IF ((IERR == 1).AND.(NPROC /= 1)) &
           & CALL ABOR1('OASIS3 grid writing only works for single CPU runs')       
        
    ENDIF
        
    ! If ierr is 1 then OASIS3 need the grid information
    ! Only write grids if we are on only one processor
        
    IF ((IERR == 1).OR.((IOASIS3LVL == 2).AND.(NPROC == 1))) THEN
        
        ALLOCATE(ZLAT(NGPTOTG,1),ZLON(NGPTOTG,1),ZLSM(NGPTOTG,1),&
          &     ZAREA(NGPTOTG,1),ILSM(NGPTOTG,1),ZCLAT(NGPTOTG,1,4),&
          &     ZCLON(NGPTOTG,1,4),ILSM2(NGPTOTG,1),ZAREA2(NGPTOTG,1))
        
        !ILSM2(:,:) = 0
        
        ! Extract the floating point LSM.
        DO JKGLO=1,NGPTOT,NPROMA
            IST =1
            IEND=MIN(NPROMA,NGPTOT-JKGLO+1)
            IBLK=(JKGLO-1)/NPROMA+1
            DO I=IST,IEND
                ZLSM(I+JKGLO-1,1)= REAL(SD_VF(I,YSD_VF%YLSM%MP,IBLK),JPRB)
            ENDDO
        ENDDO
      
        ALLOCATE (ZRGAUSLAT(0:NDGLG+1))
        ZRGAUSLAT(0)      = 90.0_JPRB
        ZRGAUSLAT(NDGLG+1)=-90.0_JPRB        
        ZRGAUSLAT(1:NDGLG)= YRCSGLEG%RLATIG(1:NDGLG)*(180.0_JPRB/(RPI))
        
        IF (LL_DEBUGCPL) THEN
           WRITE(NULOUT,'(9A12)')'No. lat','Lat.'
           DO I=0,NDGLG+1
              WRITE(NULOUT,'(I12,8F12.3)') I,ZRGAUSLAT(I)
           END DO
        END IF
        
        DO I=1,NGPTOTG
            !
            ! Setup lat,lon
            !
            ZLAT(I,1)=REAL(YRGSGEOM_NB%GELAT(I)*(180.0_JPRB/RPI),JPRB)
            ZLON(I,1)=REAL(YRGSGEOM_NB%GELAM(I)*(180.0_JPRB/RPI),JPRB)          
            !
            ! Set mask
            !
            IF (ZLSM(I,1) > RLSMCRIT) THEN
               ILSM(I,1)=1
               ILSM2(I,1)=0
            ELSE
               ILSM(I,1)=0
               ILSM2(I,1)=1
            ENDIF                        
            
            !
            ! Corners
            !
            ZDELTAX =   360.0_JPRB/REAL(NLOENG(YRGSGEOM_NB%NGPLAT(I)))
            ZDELTAYUP = 0.5_JPRB*( ZRGAUSLAT(YRGSGEOM_NB%NGPLAT(I)-1) - ZRGAUSLAT(YRGSGEOM_NB%NGPLAT(I)  ) )
            ZDELTAYDW = 0.5_JPRB*( ZRGAUSLAT(YRGSGEOM_NB%NGPLAT(I)  ) - ZRGAUSLAT(YRGSGEOM_NB%NGPLAT(I)+1) ) 
            
            ! For the northernmost latitude
            IF (ZRGAUSLAT(YRGSGEOM_NB%NGPLAT(I)-1)== 90.0_JPRB) ZDELTAYUP = 90.0_JPRB - ZRGAUSLAT(YRGSGEOM_NB%NGPLAT(I)) !2.0_JPRB*ZDELTAYUP
            ! For the southernmost latitude
            IF (ZRGAUSLAT(YRGSGEOM_NB%NGPLAT(I)+1)==-90.0_JPRB) ZDELTAYDW = ZRGAUSLAT(YRGSGEOM_NB%NGPLAT(I)) - (-90.0_JPRB) !2.0_JPRB*ZDELTAYDW
            
            !
            ! Layout of corners
            !  2 ----- 1
            !  |       |
            !  |       |
            !  3 ----- 4
            ! 
            ZCLON(I,1,1)=ZLON(I,1)+ZDELTAX*0.5_JPRB ! lon + dlon/2
            ZCLAT(I,1,1)=ZLAT(I,1)+ZDELTAYUP        ! lat + dlat/2
            ZCLON(I,1,2)=ZLON(I,1)-ZDELTAX*0.5_JPRB ! lon - dlon/2
            ZCLAT(I,1,2)=ZLAT(I,1)+ZDELTAYUP        ! lat + dlat/2
            ZCLON(I,1,3)=ZLON(I,1)-ZDELTAX*0.5_JPRB ! lon - dlon/2
            ZCLAT(I,1,3)=ZLAT(I,1)-ZDELTAYDW        ! lat - dlat/2
            ZCLON(I,1,4)=ZLON(I,1)+ZDELTAX*0.5_JPRB ! lon + dlon/2
            ZCLAT(I,1,4)=ZLAT(I,1)-ZDELTAYDW        ! lat - dlat/2
            
            WHERE (ZCLON < 0.0_JPRB) 
               ZCLON=ZCLON+360
            END WHERE
            WHERE (ZCLON > 360.0_JPRB)
               ZCLON=ZCLON-360.0_JPRB
            END WHERE           
            
            !
            ! Set areas 
            !
            ZAREA(I,1) = (ZDELTAX * RPI/180.0_JPRB * RA * COS(ZLAT(I,1)*RPI/180.0_JPRB)) * &
                       & ((ABS(ZDELTAYUP) + ABS(ZDELTAYDW)) * RPI/180.0_JPRB * RA) ![m2] 
            
            !
            ! The NEMO grid only goes down to 77.0105 degS
            ! We have to remove this area from the area calculations
            !
            
            !IF (ZCLAT(I,1,1) >= RNEMOSOUTH) THEN !if north of NEMOs southern boundary
            !   ZAREA2(I,1) = (ZDELTAX * RPI/180.0_JPRB * RA * COS(ZLAT(I,1)*RPI/180.0_JPRB)) * &
            !               & ((ABS(ZDELTAYUP) + ABS(ZDELTAYDW)) * RPI/180.0_JPRB * RA) ![m2]
            ! 
            !ELSE IF (ZCLAT(I,1,1) < RNEMOSOUTH .AND. ZCLAT(I,1,3) > RNEMOSOUTH) THEN !if southern boundary is in this latitude band
            !   ZAREA2(I,1) = (ZDELTAX * RPI/180.0_JPRB * RA * COS(ZLAT(I,1)*RPI/180.0_JPRB)) * &
            !               & (ABS(ZCLAT(I,1,1) + 77.0105) * RPI/180.0_JPRB * RA) ![m2]
            ! 
            !ELSE IF (ZCLAT(I,1,3) <= RNEMOSOUTH) THEN !if we are south of NEMOs southern boundary
            !   ZAREA2(I,1) = 0.
            !   
            !END IF            
            
            ! 
            ! Multiply grid cell area by land-sea mask
            ! to get the wet area for coupling with ocean model
            ! This includes lakes. 
            !
            ZAREA2(I,1) = ZAREA(I,1) !* (1.0_JPRB - ZLSM(I,1))
            
            IF (LL_DEBUGCPL) THEN
               IF (I==1) THEN
                  WRITE(NULOUT,'(8A12)')'No.','Lat.','Lon.','dlon','dlat', 'dx [km]','dy [km]','area [km2]'
               END IF
               WRITE(NULOUT,'(I12,7F12.3)') I,ZLAT(I,1),ZLON(I,1),ZDELTAX,ABS(ZDELTAYUP)+ABS(ZDELTAYDW),&
               & (ZDELTAX * RPI/180.0_JPRB * RA * COS(ZLAT(I,1)*RPI/180.0_JPRB))/1e3,&
               & ((ABS(ZDELTAYUP) + ABS(ZDELTAYDW)) * RPI/180.0_JPRB * RA)/1e3,ZAREA(I,1)/1e6               
            END IF 
               
        ENDDO
        
        IF (LL_DEBUGCPL) THEN
            ZGLOAREA=0.
            ZGLOAREA2=0.
            WRITE(NULOUT,'(5A9,A12,8A9)')'No.','Lat.','Lon.','LSM','ILSM', &
               &                 'Area[km2]',('Clat','Clon',J=1,4)
            DO I=1,NGPTOTG
                WRITE(NULOUT,'(I9,3F9.3,I9,F12.3,8F9.3)') &
                   & I,ZLAT(I,1),ZLON(I,1),ZLSM(I,1),ILSM(I,1),ZAREA(I,1)/1e6,&
                   & (ZCLAT(I,1,J),ZCLON(I,1,J),J=1,4)                
                ZGLOAREA = ZGLOAREA + ZAREA(I,1)            
                ZGLOAREA2 = ZGLOAREA2 + ZAREA2(I,1)
            ENDDO
            WRITE(NULOUT,'(A80,F20.3,F15.3)') 'Total global area (Google) [m2],[km2]',510.1*1e6*1e6,510.1*1e6
            WRITE(NULOUT,'(A80,F20.3,F15.3)') 'Total global area (OpenIFS)[m2],[km2]',ZGLOAREA,ZGLOAREA/1e6
            WRITE(NULOUT,'(A80,F20.3,F15.3)') 'Total global area (ORCA05)[m2],[km2]',ZGLOAREA2,ZGLOAREA2/1e6
            IF (ZGLOAREA/(510.1*1e6*1e6) > 1.02 .OR. ZGLOAREA/(510.1*1e6*1e6) < 0.98) THEN
               WRITE(NULOUT,'(A100)') 'Total global area less than 98% or more than 102% of real value'
               WRITE(NULOUT,'(A100)') 'Check how areas are calculated again '
               STOP
            ENDIF
        ENDIF
      
        ! If in true OASIS3 mode write out the data.
        IF (IOASIS3LVL == 1) THEN
            !CGRDNAM = 'atmo'
            ! Write grid
            CALL OASIS_WRITE_GRID(TRIM(CGRDNAM),NGPTOTG,1,ZLON,ZLAT)
            ! Write corners
            CALL OASIS_WRITE_CORNER(TRIM(CGRDNAM),NGPTOTG,1,4,ZCLON,ZCLAT)
            ! Write mask
            CALL OASIS_WRITE_MASK(TRIM(CGRDNAM),NGPTOTG,1,ILSM)
            ! Write areas
            CALL OASIS_WRITE_AREA(TRIM(CGRDNAM),NGPTOTG,1,ZAREA2)
            ! We are done now
            !CALL OASIS_TERMINATE_GRIDS_WRITING
            
            CGRDNAM = 'atmf'
            CALL OASIS_WRITE_GRID(TRIM(CGRDNAM),NGPTOTG,1,ZLON,ZLAT)
            ! Write corners
            CALL OASIS_WRITE_CORNER(TRIM(CGRDNAM),NGPTOTG,1,4,ZCLON,ZCLAT)
            ! Write mask
            CALL OASIS_WRITE_MASK(TRIM(CGRDNAM),NGPTOTG,1,ILSM2)
            ! Write areas
            CALL OASIS_WRITE_AREA(TRIM(CGRDNAM),NGPTOTG,1,ZAREA)
            
            CGRDNAM = 'atmr'
            ! Write grid
            CALL OASIS_WRITE_GRID(TRIM(CGRDNAM),NGPTOTG,1,ZLON,ZLAT)
            ! Write corners
            CALL OASIS_WRITE_CORNER(TRIM(CGRDNAM),NGPTOTG,1,4,ZCLON,ZCLAT)
            ! Write mask
            CALL OASIS_WRITE_MASK(TRIM(CGRDNAM),NGPTOTG,1,ILSM)
            ! Write areas
            CALL OASIS_WRITE_AREA(TRIM(CGRDNAM),NGPTOTG,1,ZAREA)
            ! We are done now
            CALL OASIS_TERMINATE_GRIDS_WRITING            
            
        ENDIF
      
        DEALLOCATE(ZLAT,ZLON,ZLSM,ZAREA,ILSM) 
        DEALLOCATE(ZRGAUSLAT)
    
    ELSE ! if we dont need to write grid files
        
        IF (LL_DEBUGCPL) THEN
            CALL CPLNG_INFO("Grid files exist, so no need to write them.")
        END IF
        
    END IF ! GWRITE == 1
    
    !END IF ! MYPROC == 1
    END IF ! LL_SKIPGRIDWRITE
    
    ! -------------------------------------------------------------------------
    ! * (6) FINALISE OASIS DEFINITION PHASE
    ! -------------------------------------------------------------------------
    WRITE(NULOUT,'(A80)') ' OpenIFS running OASIS_ENDDEF' 
    CALL OASIS_ENDDEF(IERROR)
    WRITE(NULOUT,'(A80,I8)') ' OpenIFS: OASIS_ENDDEF returns ',IERROR
    IF (IERROR/=OASIS_OK) THEN
        WRITE (CERRSTR,'(I3)') IERROR
        CALL ABOR1("CPLNG_OCE_CONFIG: Error in OASIS_ENDDEF: "//CERRSTR)
    ENDIF
    WRITE(NULOUT,'(A80)') ' OpenIFS all done in CPLNG_ADD_FLD_COMPLETED '
    
END SUBROUTINE CPLNG_ADD_FLD_COMPLETED

END MODULE CPLNG_DATA_MOD
