module RDC_core
    use fml
    use fgl
	implicit none

    !******************************************************************************
    !   data define
    !******************************************************************************
    integer(4), parameter :: nx = 128, my = 64, lev = 10
    integer(4), parameter :: LON_0 = 180
	type (mfArray), SAVE :: utd,ut,vt
	type (mfArray), SAVE :: ttd,tt,qt,pt
    real(4) :: tmax, tmin, qmax, qmin, pmax, pmin, umax, umin, vmax, vmin
    real(8) :: dat_max, dat_min

    !******************************************************************************
    !   UI flag
    !******************************************************************************
    integer(4) :: data_index, map_index, layer_index, draw_type
    logical :: animating, retry, cameraReset, updateLayer

    !******************************************************************************
    !   geo map data
    !******************************************************************************
    type (mfArray), SAVE :: x0, x1, x2, x3, y0, y1, y2, y3, z0, z1, z2, z3
    type (mfArray), SAVE :: tu0, tu1, tu2, tu3, tv0, tv1, tv2, tv3
    type (mfArray), SAVE :: ux0, ux1, ux2, ux3, uy0, uy1, uy2, uy3, uz0, uz1, uz2, uz3
    type (mfArray), SAVE :: vx0, vx1, vx2, vx3, vy0, vy1, vy2, vy3, vz0
    type (mfArray), SAVE :: cx0, cx1, cx2, cx3, cy0, cy1, cy2, cy3, cz0
    type (mfArray), SAVE :: Cam2dH,Cam3dH

contains

    subroutine CalcUnitVector( x, ux, vx )
        implicit none
        type (mfArray) :: x, y, z
        type (mfArray) :: ux, uy, uz, vx, vy, vz
        integer(4) :: dims(2)

        dims = mfShape( x )
        ux = mfS( x, MF_COL, 2.to.(dims(2)) ) - mfS( x, MF_COL, 1.to.(dims(2)-1) )
        ux = ux .hc. mfS( ux, MF_COL, 1 )
        ux = ux / ( 360.0d0 / nx ) 
        vx = mfS( x, 2.to.(dims(1)), MF_COL ) - mfS( x, 1.to.(dims(1)-1), MF_COL )
        vx = vx .vc. mfZeros( 1, dims(2) )
        vx = vx / ( 180.0d0 / (my-1) ) 
    end subroutine CalcUnitVector

    subroutine InitData()
        implicit none
        type (mfArray) :: u, v
        draw_type = 0
        map_index = 1
        data_index = 0
        layer_index = 1
        animating = .false.
        retry = .false.
		cameraReset = .true.
		updateLayer = .false.

        ! globle
        call msCreateGeoid3Data( mfOut( x0, y0, z0, tu0, tv0), LON_0-180, LON_0+180, -90, 90, nx, my-1 )
        call CalcUnitVector( x0, ux0, vx0 )
        call CalcUnitVector( y0, uy0, vy0 )
        call CalcUnitVector( z0, uz0, vz0 )

        ! transverse cylinder equal
        call msCreateGeoidData( mfOut( u, v, tu1, tv1 ), LON_0-180, LON_0+180, -90, 90, nx, my-1 )
        call msProj4( mfOut(x1, y1), 'eqc', LON_0, u, v )
        z1 = mfZeros( mfSize(x1,1),mfSize(x1,2),mfSize(x1,3))
        call CalcUnitVector( x1, ux1, vx1 )
        call CalcUnitVector( y1, uy1, vy1 )

        ! lamber Azimuth Equal Area
        call msCreateGeoidData( mfOut( u, v, tu2, tv2 ), LON_0-180, LON_0+180, 0, 90, nx, my/2-1 )
        call msProj4( mfOut(x2, y2), 'laea', LON_0, u, v, '+lat_0=90' )
        z2 = mfZeros( mfSize(x2,1),mfSize(x2,2),mfSize(x2,3))
        call CalcUnitVector( x2, ux2, vx2 )
        call CalcUnitVector( y2, uy2, vy2 )

        ! lamber Conformal Conic
        call msCreateGeoidData( mfOut( u, v, tu3, tv3 ), LON_0-180, LON_0+180, -45, 90, nx, my*3/4-1 )
        call msProj4( mfOut(x3, y3), 'lcc', LON_0, u, v, '+lat_1=20 +lat_2=60' )
        z3 = mfZeros( mfSize(x3,1),mfSize(x3,2),mfSize(x3,3))
        call CalcUnitVector( x3, ux3, vx3 )
        call CalcUnitVector( y3, uy3, vy3 )

        ! Coastline
        call msCreateCoastline3Data( mfOut( cx0, cy0, cz0 ), LON_0-180, LON_0+180, -90, 90 )
        call msCreateCoastlineData( mfOut( u, v ), LON_0-180, LON_0+180, -90, 90 )
        call msProj4( mfOut(cx1, cy1), 'eqc', LON_0, u, v )
        call msCreateCoastlineData( mfOut( u, v ), LON_0-180, LON_0+180, 0, 90 )
        call msProj4( mfOut(cx2, cy2), 'laea', LON_0, u, v, '+lat_0=90' )
        call msCreateCoastlineData( mfOut( u, v ), LON_0-180, LON_0+180, -45, 90 )
        call msProj4( mfOut(cx3, cy3), 'lcc', LON_0, u, v, '+lat_1=20 +lat_2=60' )

    end subroutine InitData
   
    subroutine LoadData()
        implicit none
		
		utd = mfLoad('masininnmi.mfb')
		ttd = mfLoad('momininnmi.mfb')
		pt = mfLoad('pt.mfb')
		
		ut = mfS(utd,MF_COL,MF_COL,MF_COL,1.to.1)
		vt = mfS(utd,MF_COL,MF_COL,MF_COL,2.to.2)
		tt = mfS(ttd,MF_COL,MF_COL,MF_COL,1.to.1)
		qt = mfS(ttd,MF_COL,MF_COL,MF_COL,2.to.2)

    end subroutine LoadData
	    
    subroutine PlotData
        if (retry) return
        retry = .true.

        call UpdateLabel

        call msHold( 'off' )

        if ( draw_type /= 2 ) then
            call PlotMasData
        end if

        if ( draw_type == 2 .or. draw_type == 3 ) then
            call PlotMomData
        end if

        if ( map_index==0 ) then
			if(cameraReset)then
				call msView('3')
				cameraReset = .false.
			end if
            call msAxis('equal')
        else
			if(cameraReset)then
				call msView('2')
				call msAxis('auto')
				call msAxis2DDependency('xy_depend',1)
				!call msAxis('tight')
				cameraReset = .false.
			end if
			if(updateLayer)then
				Cam2dH = mfGetCamViewParam('all')
				call msView('2')
				call msSetCamViewParam(mf('all'),Cam2dH)
				updateLayer=.true.
			end if
        end if
        call msDrawNow()

        retry = .false.
    end subroutine PlotData

    subroutine PlotMomData
        implicit none
        type (mfArray) :: u0, v0
        type (mfArray) :: x, y, z, tu, tv
        type (mfArray) :: u, v, w, h

		u0 = mfs(ut,MF_COL,layer_index,MF_COL)
		v0 = mfs(vt,MF_COL,layer_index,MF_COL)

        u0 = .t. mfReshape(u0,nx,my)
        v0 = .t. mfReshape(v0,nx,my)
        u0 = u0 .hc. mfS(u0,MF_COL,1)
        v0 = v0 .hc. mfS(v0,MF_COL,1)

        select case ( map_index )
        case (0)
            x = x0
            y = y0
            z = z0
            u = u0 * ux0 + v0 * vx0
            v = u0 * uy0 + v0 * vy0
            w = u0 * uz0 + v0 * vz0
            tu = tu0
            tv = tv0
        case (1)
            x = x1
            y = y1
            z = z1
            u = u0 * ux1 + v0 * vx1
            v = u0 * uy1 + v0 * vy1
            w = mfZeros( mfSize(u,1),mfSize(u,2),mfSize(u,3))
            tu = tu1
            tv = tv1
        case (2)
            x = x2
            y = y2
            z = z2
            u0 = mfS(u0, (my/2)+1 .to. my, MF_COL)
            v0 = mfS(v0, (my/2)+1 .to. my, MF_COL)
            u = u0 * ux2 + v0 * vx2
            v = u0 * uy2 + v0 * vy2
            w = mfZeros( mfSize(u,1),mfSize(u,2),mfSize(u,3))
            tu = tu2
            tv = tv2
        case (3)
            x = x3
            y = y3
            z = z3
            u0 = mfS(u0, (my/4)+1 .to. my, MF_COL)
            v0 = mfS(v0, (my/4)+1 .to. my, MF_COL)
            u = u0 * ux3 + v0 * vx3
            v = u0 * uy3 + v0 * vy3
            w = mfZeros( mfSize(u,1),mfSize(u,2),mfSize(u,3))
            tu = tu3
            tv = tv3
        end select

        if ( draw_type==2  ) then
            ! draw earth map
            h = mfSurf( x, y, z )
            call msDrawMaterial( h, 'edge', 'visible', 'off' )
            call msDrawMaterial( h, 'surf', 'smooth', 'on' )
            call msDrawTexture(h, mf('map'), mf('../data/earth.png'), &
                mf('coord_s'), tu, mf('coord_t'), tv )
	        call msDrawMaterial(h, mf('surf'), &
                mf('smooth'), mf('on'), &
		        mf('colormap'), mf('off'), &
		        mf('ambient'), mf(0), &
		        mf('diffuse'), mf(25), &
		        mf('specular'), mf(85) ) 
        end if

        call msHold( 'on' )

        if ( draw_type==2 .or. draw_type==3 ) then
            h = mfQuiver3( x, y, z, u, v, w, mf(0.1) )
            call msColorbar( 'vert' )
        end if


    end subroutine PlotMomData

    subroutine PlotMasData
        implicit none
        type (mfArray) :: x, y, z, tu, tv, cx, cy, cz
        type (mfArray) :: vdat, h

        select case ( data_index )
            case (0)   ! T
				vdat = mfS(tt,MF_COL,layer_index,MF_COL)
                dat_min = tmin
                dat_max = tmax
            case (1)   ! Q
				vdat = mfS(qt,MF_COL,layer_index,MF_COL)
            case (2)   ! P
                vdat = pt
            case (3)   ! U
				vdat = mfS(ut,MF_COL,layer_index,MF_COL)
            case (4)   ! V
				vdat = mfS(vt,MF_COL,layer_index,MF_COL)
        end select

		vdat = mfReshape(vdat,nx,my)

        vdat = .t. vdat
        vdat = vdat .hc. mfS(vdat,MF_COL,1)
        
        select case ( map_index )
        case (0)
            x = x0
            y = y0
            z = z0
            tu = tu0
            tv = tv0
            cx = cx0
            cy = cy0
            cz = cz0
        case (1)
            x = x1
            y = y1
            z = z1
            tu = tu1
            tv = tv1
            cx = cx1
            cy = cy1
        case (2)
            vdat = mfS(vdat, (my/2)+1 .to. my, MF_COL)
            x = x2
            y = y2
            z = z2
            tu = tu2
            tv = tv2
            cx = cx2
            cy = cy2
        case (3)
            vdat = mfS(vdat, (my/4)+1 .to. my, MF_COL)
            x = x3
            y = y3
            z = z3
            tu = tu3
            tv = tv3
            cx = cx3
            cy = cy3
        end select
       
        dat_max = mfMax( mfS( vdat, MF_COL ) )
        dat_min = mfMin( mfS( vdat, MF_COL ) )
            
        if ( draw_type==0 .or. draw_type==3 ) then
            ! draw earth map
            h = mfSurf( x, y, z )
            call msDrawMaterial( h, 'edge', 'visible', 'off' )
            call msDrawMaterial( h, 'surf', 'smooth', 'on' )
            call msDrawTexture(h, mf('map'), mf('../data/earth.png'), &
                mf('coord_s'), tu, mf('coord_t'), tv )
	        call msDrawMaterial(h, mf('surf'), &
                mf('smooth'), mf('on'), &
		        mf('colormap'), mf('off'), &
		        mf('ambient'), mf(0), &
		        mf('diffuse'), mf(25), &
		        mf('specular'), mf(85) ) 
            call msHold( 'on' )
        end if
        
        if ( draw_type==0 .or. draw_type==1 ) then
            if ( map_index==0 ) then
                h = mfPlot3( cx, cy, cz )
                call msDrawMaterial( h, 'edge', 'colormap', 'off' )
            else
                h = mfPlot( cx, cy )
            end if
            call msHold( 'on' )
        end if


        if ( draw_type==0 .or. draw_type==3 ) then
            h = mfContour3( x, y, z, vdat )
        else if ( draw_type==1 ) then           
            h = mfSolidContour3( x, y, z, vdat )
            call msDrawMaterial( h, 'edge', 'visible', 'off' )
            call msDrawMaterial( h, 'surf', 'smooth', 'on' )
        end if
        call msColorbar( 'vert' )
        call msColormapRange( dat_min, dat_max )  
    end subroutine PlotMasData

    subroutine UpdateLabel
        use mxui
        implicit none
	    character(80) :: ch

	    write (ch,100) layer_index
        call msUISetPropertyString( 'lblLayer', 'caption', ch ) 
        call msUISetPropertyInteger( 'scbLayer', 'position', layer_index ) 

    100 format("Layer: (",i2,"/10)")

    end subroutine UpdateLabel

    subroutine Animation

        do while (animating) 
            if (layer_index>=lev) exit
            layer_index = layer_index + 1
            call PlotData
        end do

        animating = .false.

    end subroutine Animation

    subroutine StartAnimation
        if ( animating ) return
        if ( data_index==2 ) return  ! P
        animating = .true.
        layer_index = 0
        call Animation
    end subroutine StartAnimation

    subroutine StopAnimation
        animating = .false.
    end subroutine StopAnimation


end module RDC_core