MODULE etat0_database_mod


CONTAINS

  SUBROUTINE init_etat0
  USE xios_mod
  USE omp_para
  IMPLICIT NONE
  
    IF (is_omp_master) THEN
      CALL xios_set_fieldgroup_attr("read_fields",enabled=.TRUE.)
      CALL xios_set_filegroup_attr("read_files",enabled=.TRUE.)
    ENDIF
  END SUBROUTINE init_etat0

  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
  USE icosa
  USE restart_mod
  USE wind_mod
  USE write_field_mod
  USE time_mod
  USE transfert_mod
  USE xios_mod
  USE write_field_mod
  USE vertical_remap_mod
  USE theta2theta_rhodz_mod
  USE qsat_mod
  USE pression_mod
  USE omp_para
  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_ulon_reg(:)
    TYPE(t_field),POINTER,SAVE :: f_ulat_reg(:)
    TYPE(t_field),POINTER,SAVE :: f_temp_reg(:)
    TYPE(t_field),POINTER,SAVE :: f_q_reg(:)

    TYPE(t_field),POINTER,SAVE :: f_ts(:)
    TYPE(t_field),POINTER,SAVE :: f_z(:)
    TYPE(t_field),POINTER,SAVE :: f_ulon(:)
    TYPE(t_field),POINTER,SAVE :: f_ulat(:)
    TYPE(t_field),POINTER,SAVE :: f_temp(:)
    TYPE(t_field),POINTER,SAVE :: f_q1(:)
    TYPE(t_field),POINTER,SAVE :: f_qsat(:)
    TYPE(t_field),POINTER,SAVE :: f_p(:)
    INTEGER :: nb_level
    REAL,ALLOCATABLE:: levels(:)
    INTEGER :: ind

    CALL xios_read_field("relief_db",f_phis)
    
    CALL writeField("relief_out",f_phis,once=.TRUE.)

    IF (is_omp_level_master) THEN
      DO ind=1,ndomain
        IF (.NOT. assigned_domain(ind)) CYCLE
        f_phis(ind)%rval2d(:)=f_phis(ind)%rval2d(:)*g      
      ENDDO
    ENDIF
!$OMP BARRIER    

    IF (is_omp_master) CALL xios_get_axis_attr("lev_ecdyn",n_glo=nb_level)
    CALL bcast_omp(nb_level)
    ALLOCATE(levels(nb_level))

    IF (is_omp_master) CALL xios_get_axis_attr("lev_ecdyn",value=levels)
    CALL bcast_omp(levels)
    
    levels=levels*100  ! hectoPascal -> Pascal
 
    CALL allocate_field(f_ts, field_t, type_real, name="ts")
    CALL allocate_field(f_z, field_t, type_real, name="z")
    CALL allocate_field(f_ulon_reg, field_t, type_real,nb_level)
    CALL allocate_field(f_ulat_reg, field_t, type_real,nb_level)
    CALL allocate_field(f_temp_reg, field_t, type_real,nb_level)
    CALL allocate_field(f_q_reg,    field_t, type_real,nb_level)
    
    CALL allocate_field(f_q1,    field_t, type_real,llm)
    CALL allocate_field(f_qsat,  field_t, type_real,llm)
    CALL allocate_field(f_p,  field_t, type_real,llm+1)
    CALL allocate_field(f_temp,  field_t, type_real,llm)
    CALL allocate_field(f_ulon,  field_t, type_real,llm)
    CALL allocate_field(f_ulat,  field_t, type_real,llm)

    CALL xios_read_field("z_db",f_z)
    CALL xios_read_field("ps_db",f_ps)
    CALL xios_read_field("ts_db",f_ts)
    CALL writeField("ps_out",f_ps)

!$OMP BARRIER
    
!    CALL writeField("phis_out",f_phis,once=.TRUE.)
!    CALL writeField("ts_out",f_ts,once=.TRUE.)

! make correction to ps due to relief at higher resolution
! difference with LMDZ : tsol is taken from ECDYN.NC and not from ECPHY
    IF (is_omp_level_master) THEN
      DO ind=1,ndomain
        IF (.NOT. assigned_domain(ind)) CYCLE
        f_ps(ind)%rval2d(:)=f_ps(ind)%rval2d(:)*(1.+(f_z(ind)%rval2d(:)-f_phis(ind)%rval2d(:))/287.0/f_ts(ind)%rval2d(:))
      ENDDO
    ENDIF
!$OMP BARRIER    
    CALL transfert_request(f_ps,req_i0) 
    CALL writeField("ps_out",f_ps)
    
    

    CALL xios_read_field("temp_db",f_temp_reg)
    CALL vertical_remap(levels,f_temp_reg,f_ps,f_temp)
    CALL transfert_request(f_temp,req_i0) 
    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)

    CALL xios_read_field("u_db",f_ulon_reg)
    CALL vertical_remap(levels,f_ulon_reg,f_ps,f_ulon)
    CALL xios_read_field("v_db",f_ulat_reg)
    CALL vertical_remap(levels,f_ulat_reg,f_ps,f_ulat)
    CALL transfert_request(f_ulat,req_i0) 
    CALL transfert_request(f_ulon,req_i0) 
    CALL ulonlat2un(f_ulon, f_ulat,f_u)

    CALL xios_read_field("q_db",f_q_reg)
    CALL vertical_remap(levels,f_q_reg,f_ps,f_q1)

    CALL pression(f_ps,f_p)
! difference with LMDZ : for qsat, pressure at mid layer is computed as a mean value pmid=(p(l)+p(l+1))/2    
    CALL qsat(f_temp,f_p,f_qsat)
    CALL transfert_request(f_qsat,req_i0) 

    DO ind=1,ndomain
      IF (.NOT. assigned_domain(ind)) CYCLE
      f_q(ind)%rval4d(:,:,:)=1e-6
      f_q(ind)%rval4d(:,:,1)=f_q1(ind)%rval3d(:,:)*f_qsat(ind)%rval3d(:,:)*0.01
      WHERE(f_q(ind)%rval4d(:,:,1)<0) f_q(ind)%rval4d(:,:,1)=0
    ENDDO

    CALL writeField("tempdb_out",f_temp_reg)
    CALL writeField("temp_out",f_temp)

    CALL deallocate_field(f_ts)
    CALL deallocate_field(f_z)
    CALL deallocate_field(f_ulon_reg)
    CALL deallocate_field(f_ulat_reg)
    CALL deallocate_field(f_temp_reg)
    CALL deallocate_field(f_q_reg)
    
    CALL deallocate_field(f_q1)
    CALL deallocate_field(f_qsat)
    CALL deallocate_field(f_p)
    CALL deallocate_field(f_temp)
    CALL deallocate_field(f_ulon)
    CALL deallocate_field(f_ulat)   
 
  END SUBROUTINE etat0

END MODULE etat0_database_mod