MODULE etat0_dcmip2_mod

! test cases DCMIP 2012, category 2 : Orographic gravity waves

  USE genmod

  PRIVATE

  PUBLIC  etat0

CONTAINS


  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
  USE icosa
  USE theta2theta_rhodz_mod
  IMPLICIT NONE
    TYPE(t_field),POINTER :: f_ps(:)
    TYPE(t_field),POINTER :: f_phis(:)
    TYPE(t_field),POINTER :: f_theta_rhodz(:)
    TYPE(t_field),POINTER :: f_u(:)
    TYPE(t_field),POINTER :: f_q(:)
        
    TYPE(t_field),POINTER,SAVE :: f_temp(:)
    REAL(rstd),POINTER :: ps(:)
    REAL(rstd),POINTER :: phis(:)
    REAL(rstd),POINTER :: u(:,:)
    REAL(rstd),POINTER :: Temp(:,:)
    REAL(rstd),POINTER :: q(:,:,:)

    INTEGER :: ind, icase
    CHARACTER(len=255) :: etat0_type
    etat0_type='jablonowsky06'
    CALL getin("etat0",etat0_type)

    SELECT CASE (TRIM(etat0_type))
    CASE('dcmip2_mountain')
       icase=0
    CASE('dcmip2_schaer_noshear')
       icase=1
    CASE('dcmip2_schaer_shear')
       icase=2
    CASE DEFAULT
       PRINT *, 'This should not happen : etat0_type =', TRIM(etat0_type), ' in etat0_dcmip2.f90/etat0'
       STOP
    END SELECT
    PRINT *, 'Orographic gravity-wave test case :', TRIM(etat0_type)

    CALL allocate_field(f_temp,field_t,type_real,llm,name='temp') 

    DO ind=1,ndomain
      IF (.NOT. assigned_domain(ind)) CYCLE
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)

      ps=f_ps(ind)
      phis=f_phis(ind)
      u=f_u(ind)
      q=f_q(ind)
      temp=f_temp(ind)
      CALL compute_etat0_DCMIP2(icase,ps,phis,u,Temp,q)
    ENDDO
    
    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
            
    CALL deallocate_field(f_temp)
   
  END SUBROUTINE etat0
  
  SUBROUTINE compute_etat0_DCMIP2(icase, ps, phis, u, temp,q)
    USE icosa
    USE disvert_mod, ONLY : ap,bp
    USE pression_mod
    USE theta2theta_rhodz_mod
    USE wind_mod
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: icase
    REAL(rstd), INTENT(OUT) :: ps(iim*jjm)
    REAL(rstd), INTENT(OUT) :: phis(iim*jjm)
    REAL(rstd), INTENT(OUT) :: u(3*iim*jjm,llm)
    REAL(rstd), INTENT(OUT) :: temp(iim*jjm,llm)
    REAL(rstd), INTENT(OUT) :: q(iim*jjm,llm)
    REAL(rstd) :: ulon(3*iim*jjm,llm), ulat(3*iim*jjm,llm)
    
    INTEGER :: i,j,l,ij
    REAL(rstd) :: hyam, hybm, pp, dummy1, dummy2, dummy3  
    
    ! Hexagons : ps,phis,temp
    DO l=1,llm
       ! The surface pressure is not set yet so we provide the hybrid coefficients
       hyam = .5*(ap(l)+ap(l+1))/preff
       hybm = .5*(bp(l)+bp(l+1))
       DO j=jj_begin,jj_end
          DO i=ii_begin,ii_end
             ij=(j-1)*iim+i
             CALL comp_all(lon_i(ij),lat_i(ij), ps(ij),phis(ij),temp(ij,l), dummy1,dummy2)
          END DO
       END DO
    END DO
    
    ! Edges : ulon,ulat
    DO l=1,llm
       ! The surface pressure is not set yet so we provide the hybrid coefficients
       hyam = .5*(ap(l)+ap(l+1))/preff
       hybm = .5*(bp(l)+bp(l+1))
       DO j=jj_begin,jj_end
          DO i=ii_begin,ii_end
             ij=(j-1)*iim+i
             CALL comp_all(lon_e(ij+u_right), lat_e(ij+u_right), &
                  dummy1,dummy2,dummy3, ulon(ij+u_right,l),ulat(ij+u_right,l))
             CALL comp_all(lon_e(ij+u_lup), lat_e(ij+u_lup), &
                  dummy1,dummy2,dummy3, ulon(ij+u_lup,l),ulat(ij+u_lup,l))
             CALL comp_all(lon_e(ij+u_ldown), lat_e(ij+u_ldown), &
                  dummy1,dummy2,dummy3, ulon(ij+u_ldown,l),ulat(ij+u_ldown,l))
          END DO
       END DO
    END DO
    CALL compute_wind_perp_from_lonlat_compound(ulon,ulat,u)

    q=1.

    CONTAINS
      
      SUBROUTINE comp_all(lon,lat, psj,phisj,tempj, ulonj,ulatj)
        USE dcmip_initial_conditions_test_1_2_3
        REAL(rstd), INTENT(IN) :: lon, lat
        REAL(rstd), INTENT(OUT) :: psj,phisj,tempj,ulonj,ulatj
        REAL :: dummy
        dummy=0.
        SELECT CASE (icase)
        CASE(0)
           CALL test2_steady_state_mountain(lon,lat,dummy,dummy,0,.TRUE.,hyam,hybm, &
                ulonj,ulatj,dummy,tempj,phisj,psj,dummy,dummy)
        CASE(1)
           CALL test2_schaer_mountain(lon,lat,dummy,dummy,0,.TRUE.,hyam,hybm,0,&
                ulonj,ulatj,dummy,tempj,phisj,psj,dummy,dummy)
        CASE(2)
           CALL test2_schaer_mountain(lon,lat,dummy,dummy,0,.TRUE.,hyam,hybm,1, &
                ulonj,ulatj,dummy,tempj,phisj,psj,dummy,dummy)
        END SELECT
      END SUBROUTINE comp_all

  END SUBROUTINE compute_etat0_DCMIP2


END MODULE etat0_dcmip2_mod