MODULE lmdz_call_calfis

  USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_DEBUGIO

  REAL, POINTER, SAVE :: ucov(:, :)
  REAL, POINTER, SAVE :: vcov(:, :)
  REAL, POINTER, SAVE :: teta(:, :)
  REAL, POINTER, SAVE :: masse(:, :)
  REAL, POINTER, SAVE :: ps(:)
  REAL, POINTER, SAVE :: phis(:)
  REAL, POINTER, SAVE :: q(:, :, :)
  REAL, POINTER, SAVE :: flxw(:, :)

  REAL, POINTER, SAVE :: p(:, :)
  REAL, POINTER, SAVE :: pks(:)
  REAL, POINTER, SAVE :: pk(:, :)
  REAL, POINTER, SAVE :: pkf(:, :)
  REAL, POINTER, SAVE :: phi(:, :)
  REAL, POINTER, SAVE :: du(:, :)
  REAL, POINTER, SAVE :: dv(:, :)
  REAL, POINTER, SAVE :: dteta(:, :)
  REAL, POINTER, SAVE :: dq(:, :, :)
  REAL, POINTER, SAVE :: dufi(:, :)
  REAL, POINTER, SAVE :: dvfi(:, :)
  REAL, POINTER, SAVE :: dtetafi(:, :)
  REAL, POINTER, SAVE :: dqfi(:, :, :)
  REAL, POINTER, SAVE :: dpfi(:)


CONTAINS

  SUBROUTINE call_calfis_allocate
    USE bands
    USE allocate_field_mod
    USE parallel_lmdz
    USE lmdz_dimensions
    USE lmdz_paramet
    USE infotrac, ONLY: nqtot
    IMPLICIT NONE
    TYPE(distrib), POINTER :: d
    d => distrib_physic

    CALL allocate_u(ucov, llm, d)
    CALL allocate_v(vcov, llm, d)
    CALL allocate_u(teta, llm, d)
    CALL allocate_u(masse, llm, d)
    CALL allocate_u(ps, d)
    CALL allocate_u(phis, d)
    CALL allocate_u(q, llm, nqtot, d)
    CALL allocate_u(flxw, llm, d)
    CALL allocate_u(p, llmp1, d)
    CALL allocate_u(pks, d)
    pks(:) = 0
    CALL allocate_u(pk, llm, d)
    pk(:, :) = 0
    CALL allocate_u(pkf, llm, d)
    CALL allocate_u(phi, llm, d)
    CALL allocate_u(du, llm, d)
    CALL allocate_v(dv, llm, d)
    CALL allocate_u(dteta, llm, d)
    CALL allocate_u(dq, llm, nqtot, d)
    CALL allocate_u(dufi, llm, d)
    CALL allocate_v(dvfi, llm, d)
    CALL allocate_u(dtetafi, llm, d)
    CALL allocate_u(dqfi, llm, nqtot, d)
    CALL allocate_u(dpfi, d)

  END SUBROUTINE call_calfis_allocate


  SUBROUTINE call_calfis(itau, lafin, ucov_dyn, vcov_dyn, teta_dyn, masse_dyn, ps_dyn, &
          phis_dyn, q_dyn, flxw_dyn)
    USE lmdz_dimensions
    USE lmdz_paramet
    USE exner_hyb_loc_m, ONLY: exner_hyb_loc
    USE exner_milieu_loc_m, ONLY: exner_milieu_loc
    USE parallel_lmdz
    USE times
    USE mod_hallo
    USE Bands
    USE lmdz_vampir
    USE infotrac, ONLY: nqtot
    USE control_mod
    USE write_field_loc
    USE lmdz_strings, ONLY: int2str
    USE comconst_mod, ONLY: dtphys
    USE logic_mod, ONLY: leapf, forward, ok_strato
    USE comvert_mod, ONLY: ap, bp, pressure_exner
    USE temps_mod, ONLY: day_ini, day_ref, jd_ref, jh_ref, start_time
    USE lmdz_calfis_loc
    USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS
    USE lmdz_iniprint, ONLY: lunout, prt_level

    IMPLICIT NONE

    INTEGER, INTENT(IN) :: itau ! (time) iteration step number
    LOGICAL, INTENT(IN) :: lafin ! .TRUE. if final time step
    REAL, INTENT(INOUT) :: ucov_dyn(ijb_u:ije_u, llm) ! covariant zonal wind
    REAL, INTENT(INOUT) :: vcov_dyn(ijb_v:ije_v, llm) ! covariant meridional wind
    REAL, INTENT(INOUT) :: teta_dyn(ijb_u:ije_u, llm) ! potential temperature
    REAL, INTENT(INOUT) :: masse_dyn(ijb_u:ije_u, llm) ! air mass
    REAL, INTENT(INOUT) :: ps_dyn(ijb_u:ije_u) ! surface pressure
    REAL, INTENT(INOUT) :: phis_dyn(ijb_u:ije_u) ! surface geopotential
    REAL, INTENT(INOUT) :: q_dyn(ijb_u:ije_u, llm, nqtot) ! advected tracers
    REAL, INTENT(INOUT) :: flxw_dyn(ijb_u:ije_u, llm) ! vertical mass flux

    REAL :: dufi_tmp(iip1, llm)
    REAL :: dvfi_tmp(iip1, llm)
    REAL :: dtetafi_tmp(iip1, llm)
    REAL :: dpfi_tmp(iip1)
    REAL :: dqfi_tmp(iip1, llm, nqtot)

    REAL :: jD_cur, jH_cur
    CHARACTER(LEN = 15) :: ztit
    TYPE(Request), SAVE :: Request_physic
    !$OMP THREADPRIVATE(Request_physic )
    INTEGER :: ijb, ije, l, iq

    IF (CPPKEY_DEBUGIO) THEN
      CALL WriteField_u('ucovfi', ucov)
      CALL WriteField_v('vcovfi', vcov)
      CALL WriteField_u('tetafi', teta)
      CALL WriteField_u('pfi', p)
      CALL WriteField_u('pkfi', pk)
      DO iq = 1, nqtot
        CALL WriteField_u('qfi' // trim(int2str(iq)), q(:, :, iq))
      ENDDO
    END IF

    !     .......   Ajout   P.Le Van ( 17/04/96 )   ...........

    !$OMP MASTER
    CALL suspend_timer(timer_caldyn)
    IF (prt_level >= 10) THEN
      WRITE(lunout, *) 'leapfrog_p: Entree dans la physique : Iteration No ', itau
    ENDIF
    !$OMP END MASTER

    jD_cur = jD_ref + day_ini - day_ref + (itau + 1) / day_step

    IF (planet_type =="generic") THEN
      ! AS: we make jD_cur to be pday
      jD_cur = int(day_ini + itau / day_step)
    ENDIF

    jH_cur = jH_ref + start_time + &
            mod(itau + 1, day_step) / float(day_step)
    IF (jH_cur > 1.0) THEN
      jD_cur = jD_cur + 1.
      jH_cur = jH_cur - 1.
    endif

    !   Inbterface avec les routines de phylmd (phymars ... )
    !   -----------------------------------------------------

    !+jld

    !  Diagnostique de conservation de l'energie : initialisation

    !-jld
    !$OMP BARRIER
    !$OMP MASTER
    CALL VTb(VThallo)
    !$OMP END MASTER

    IF (CPPKEY_DEBUGIO) THEN
      CALL WriteField_u('ucovfi', ucov)
      CALL WriteField_v('vcovfi', vcov)
      CALL WriteField_u('tetafi', teta)
      CALL WriteField_u('pfi', p)
      CALL WriteField_u('pkfi', pk)
    END IF

    CALL SetTag(Request_physic, 800)
    CALL Register_SwapField_u(ucov_dyn, ucov, distrib_physic, Request_physic, up = 2, down = 2)
    CALL Register_SwapField_v(vcov_dyn, vcov, distrib_physic, Request_physic, up = 2, down = 2)
    CALL Register_SwapField_u(teta_dyn, teta, distrib_physic, Request_physic, up = 2, down = 2)
    CALL Register_SwapField_u(masse_dyn, masse, distrib_physic, Request_physic, up = 1, down = 2)
    CALL Register_SwapField_u(ps_dyn, ps, distrib_physic, Request_physic, up = 2, down = 2)
    CALL Register_SwapField_u(phis_dyn, phis, distrib_physic, Request_physic, up = 2, down = 2)
    CALL Register_SwapField_u(q_dyn, q, distrib_physic, Request_physic, up = 2, down = 2)
    CALL Register_SwapField_u(flxw_dyn, flxw, distrib_physic, Request_physic, up = 2, down = 2)

    CALL SendRequest(Request_Physic)
    !$OMP BARRIER
    CALL WaitRequest(Request_Physic)

    !$OMP BARRIER
    !$OMP MASTER
    CALL Set_Distrib(distrib_Physic)
    CALL VTe(VThallo)

    CALL VTb(VTphysiq)
    !$OMP END MASTER
    !$OMP BARRIER

    CALL pression_loc (ip1jmp1, ap, bp, ps, p)

    !$OMP BARRIER
    CALL exner_hyb_loc(ip1jmp1, ps, p, pks, pk, pkf)
    !$OMP BARRIER
    CALL geopot_loc  (ip1jmp1, teta, pk, pks, phis, phi)

    CALL Register_Hallo_u(p, llmp1, 2, 2, 2, 2, Request_physic)
    CALL Register_Hallo_u(pk, llm, 2, 2, 2, 2, Request_physic)
    CALL Register_Hallo_u(phi, llm, 2, 2, 2, 2, Request_physic)

    CALL SendRequest(Request_Physic)
    !$OMP BARRIER
    CALL WaitRequest(Request_Physic)

    !$OMP BARRIER

    IF (CPPKEY_DEBUGIO) THEN
      CALL WriteField_u('ucovfi', ucov)
      CALL WriteField_v('vcovfi', vcov)
      CALL WriteField_u('tetafi', teta)
      CALL WriteField_u('pfi', p)
      CALL WriteField_u('pkfi', pk)
      DO iq = 1, nqtot
        CALL WriteField_u('qfi' // trim(int2str(iq)), q(:, :, iq))
      ENDDO
    END IF

    !$OMP BARRIER

    IF (CPPKEY_PHYS) THEN
      CALL calfis_loc(lafin, jD_cur, jH_cur, &
              ucov, vcov, teta, q, masse, ps, p, pk, phis, phi, &
              du, dv, dteta, dq, &
              flxw, dufi, dvfi, dtetafi, dqfi, dpfi)
    END IF
    ijb = ij_begin
    ije = ij_end
    IF (.NOT. pole_nord) THEN

      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l = 1, llm
        dufi_tmp(1:iip1, l) = dufi(ijb:ijb + iim, l)
        dvfi_tmp(1:iip1, l) = dvfi(ijb:ijb + iim, l)
        dtetafi_tmp(1:iip1, l) = dtetafi(ijb:ijb + iim, l)
        dqfi_tmp(1:iip1, l, :) = dqfi(ijb:ijb + iim, l, :)
      ENDDO
      !$OMP END DO NOWAIT

      !$OMP MASTER
      dpfi_tmp(1:iip1) = dpfi(ijb:ijb + iim)
      !$OMP END MASTER

    ENDIF ! of if ( .NOT. pole_nord)

    !$OMP BARRIER
    !$OMP MASTER
    CALL Set_Distrib(distrib_Physic_bis)
    CALL VTb(VThallo)
    !$OMP END MASTER
    !$OMP BARRIER

    CALL Register_Hallo_u(dufi, llm, 1, 0, 0, 1, Request_physic)
    CALL Register_Hallo_v(dvfi, llm, 1, 0, 0, 1, Request_physic)
    CALL Register_Hallo_u(dtetafi, llm, 1, 0, 0, 1, Request_physic)
    CALL Register_Hallo_u(dpfi, 1, 1, 0, 0, 1, Request_physic)

    DO iq = 1, nqtot
      CALL Register_Hallo_u(dqfi(:, :, iq), llm, 1, 0, 0, 1, Request_physic)
    ENDDO

    CALL SendRequest(Request_Physic)
    !$OMP BARRIER
    CALL WaitRequest(Request_Physic)

    !$OMP BARRIER
    !$OMP MASTER
    CALL VTe(VThallo)
    CALL Set_Distrib(distrib_Physic)
    !$OMP END MASTER
    !$OMP BARRIER
    ijb = ij_begin
    IF (.NOT. pole_nord) THEN

      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l = 1, llm
        dufi(ijb:ijb + iim, l) = dufi(ijb:ijb + iim, l) + dufi_tmp(1:iip1, l)
        dvfi(ijb:ijb + iim, l) = dvfi(ijb:ijb + iim, l) + dvfi_tmp(1:iip1, l)
        dtetafi(ijb:ijb + iim, l) = dtetafi(ijb:ijb + iim, l) + dtetafi_tmp(1:iip1, l)
        dqfi(ijb:ijb + iim, l, :) = dqfi(ijb:ijb + iim, l, :) + dqfi_tmp(1:iip1, l, :)
      ENDDO
      !$OMP END DO NOWAIT

      !$OMP MASTER
      dpfi(ijb:ijb + iim) = dpfi(ijb:ijb + iim) + dpfi_tmp(1:iip1)
      !$OMP END MASTER

    endif ! of if (.NOT. pole_nord)

    IF (CPPKEY_DEBUGIO) THEN
      CALL WriteField_u('dufi', dufi)
      CALL WriteField_v('dvfi', dvfi)
      CALL WriteField_u('dtetafi', dtetafi)
      CALL WriteField_u('dpfi', dpfi)
      DO iq = 1, nqtot
        CALL WriteField_u('dqfi' // trim(int2str(iq)), dqfi(:, :, iq))
      ENDDO
    END IF

    !$OMP BARRIER

    !      ajout des tendances physiques:
    !      ------------------------------
    IF (CPPKEY_DEBUGIO) THEN
      CALL WriteField_u('ucovfi', ucov)
      CALL WriteField_v('vcovfi', vcov)
      CALL WriteField_u('tetafi', teta)
      CALL WriteField_u('psfi', ps)
      DO iq = 1, nqtot
        CALL WriteField_u('qfi' // trim(int2str(iq)), q(:, :, iq))
      ENDDO
    END IF

    IF (CPPKEY_DEBUGIO) THEN
      CALL WriteField_u('ucovfi', ucov)
      CALL WriteField_v('vcovfi', vcov)
      CALL WriteField_u('tetafi', teta)
      CALL WriteField_u('psfi', ps)
      DO iq = 1, nqtot
        CALL WriteField_u('qfi' // trim(int2str(iq)), q(:, :, iq))
      ENDDO
    END IF

    CALL addfi_loc(dtphys, leapf, forward, &
            ucov, vcov, teta, q, ps, &
            dufi, dvfi, dtetafi, dqfi, dpfi)
    ! since addfi updates ps(), also update p(), masse() and pk()
    CALL pression_loc(ip1jmp1, ap, bp, ps, p)
    !$OMP BARRIER
    CALL massdair_loc(p, masse)
    !$OMP BARRIER
    IF (pressure_exner) THEN
      CALL exner_hyb_loc(ijnb_u, ps, p, pks, pk, pkf)
    else
      CALL exner_milieu_loc(ijnb_u, ps, p, pks, pk, pkf)
    endif
    !$OMP BARRIER

    IF (CPPKEY_DEBUGIO) THEN
      CALL WriteField_u('ucovfi', ucov)
      CALL WriteField_v('vcovfi', vcov)
      CALL WriteField_u('tetafi', teta)
      CALL WriteField_u('psfi', ps)
      DO iq = 1, nqtot
        CALL WriteField_u('qfi' // trim(int2str(iq)), q(:, :, iq))
      ENDDO
    END IF

    IF (ok_strato) THEN
      !      CALL top_bound_loc( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
      CALL top_bound_loc(vcov, ucov, teta, masse, dtphys)
    ENDIF

    !$OMP BARRIER
    !$OMP MASTER
    CALL VTe(VTphysiq)
    CALL VTb(VThallo)
    !$OMP END MASTER

    CALL SetTag(Request_physic, 800)
    CALL Register_SwapField_u(ucov, ucov_dyn, distrib_caldyn, Request_physic)
    CALL Register_SwapField_v(vcov, vcov_dyn, distrib_caldyn, Request_physic)
    CALL Register_SwapField_u(teta, teta_dyn, distrib_caldyn, Request_physic)
    CALL Register_SwapField_u(masse, masse_dyn, distrib_caldyn, Request_physic)
    CALL Register_SwapField_u(ps, ps_dyn, distrib_caldyn, Request_physic)
    CALL Register_SwapField_u(q, q_dyn, distrib_caldyn, Request_physic)
    CALL SendRequest(Request_Physic)
    !$OMP BARRIER
    CALL WaitRequest(Request_Physic)

    !$OMP BARRIER
    !$OMP MASTER
    CALL VTe(VThallo)
    CALL set_distrib(distrib_caldyn)
    !$OMP END MASTER
    !$OMP BARRIER

    !  Diagnostique de conservation de l'energie : difference
    IF (ip_ebil_dyn>=1) THEN
      ztit = 'bil phys'
      !      CALL diagedyn(ztit,2,1,1,dtphys,ucov, vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
      WRITE(lunout, *)"call_calfis: diagedyn disabled in dyn3dmem !!"
    ENDIF

    IF (CPPKEY_DEBUGIO) THEN
      CALL WriteField_u('ucovfi', ucov_dyn)
      CALL WriteField_v('vcovfi', vcov_dyn)
      CALL WriteField_u('tetafi', teta_dyn)
      CALL WriteField_u('psfi', ps_dyn)
      DO iq = 1, nqtot
        CALL WriteField_u('qfi' // trim(int2str(iq)), q_dyn(:, :, iq))
      ENDDO
    END IF


    !-jld
    !$OMP MASTER
    CALL resume_timer(timer_caldyn)
    !$OMP END MASTER

  END SUBROUTINE call_calfis

END MODULE lmdz_call_calfis