source: LMDZ6/trunk/libf/phylmd/create_etat0_unstruct_mod.F90 @ 5134

MODULE create_etat0_unstruct_mod

  REAL, SAVE, ALLOCATABLE :: zmea_gw(:)
  !$OMP THREADPRIVATE(zmea_gw)
  REAL, SAVE, ALLOCATABLE :: zpic_gw(:)
  !$OMP THREADPRIVATE(zpic_gw)
  REAL, SAVE, ALLOCATABLE :: zval_gw(:)
  !$OMP THREADPRIVATE(zval_gw)
  REAL, SAVE, ALLOCATABLE :: zstd_gw(:)
  !$OMP THREADPRIVATE(zstd_gw)
  REAL, SAVE, ALLOCATABLE :: zsig_gw(:)
  !$OMP THREADPRIVATE(zsig_gw)
  REAL, SAVE, ALLOCATABLE :: zgam_gw(:)
  !$OMP THREADPRIVATE(zgam_gw)
  REAL, SAVE, ALLOCATABLE :: zthe_gw(:)
  !$OMP THREADPRIVATE(zthe_gw)

  PRIVATE zmea_gw, zpic_gw, zval_gw, zstd_gw, zsig_gw, zgam_gw, zthe_gw


CONTAINS
  
  SUBROUTINE init_create_etat0_unstruct
  USE lmdz_xios
  USE netcdf
  USE mod_phys_lmdz_para
  IMPLICIT NONE
  INTEGER :: file_id, iret
  
   ! for coupling activate ocean fraction reading from file "ocean_fraction.nc"
    IF (is_omp_master) THEN

      IF (NF90_OPEN("ocean_fraction.nc", NF90_NOWRITE, file_id)==NF90_NOERR) THEN
        CALL xios_set_file_attr("frac_ocean",enabled=.TRUE.)
        CALL xios_set_field_attr("mask",field_ref="frac_ocean_read")
        iret=NF90_CLOSE(file_id)
      ELSE IF (NF90_OPEN("land_water_0.05.nc", NF90_NOWRITE, file_id)==NF90_NOERR) THEN
        CALL xios_set_file_attr("land_water",name="land_water_0.05",enabled=.TRUE.)
        CALL xios_set_field_attr("mask",field_ref="land_water")
        iret=NF90_CLOSE(file_id)
      ELSE IF (NF90_OPEN("land_water_0.25.nc", NF90_NOWRITE, file_id)==NF90_NOERR) THEN
        CALL xios_set_file_attr("land_water",name="land_water_0.25",enabled=.TRUE.)
        CALL xios_set_field_attr("mask",field_ref="land_water")
        iret=NF90_CLOSE(file_id)
      ELSE IF (NF90_OPEN("land_water_0.50.nc", NF90_NOWRITE, file_id)==NF90_NOERR) THEN
        CALL xios_set_file_attr("land_water",name="land_water_0.50",enabled=.TRUE.)
        CALL xios_set_field_attr("mask",field_ref="land_water")
        iret=NF90_CLOSE(file_id)
      ENDIF

    ENDIF

  END SUBROUTINE init_create_etat0_unstruct
  
 
  SUBROUTINE init_param_gw(zmea, zpic, zval, zstd, zsig, zgam, zthe)
  USE dimphy
    REAL, INTENT(IN) :: zmea(klon)  
    REAL, INTENT(IN) :: zpic(klon)
    REAL, INTENT(IN) :: zval(klon) 
    REAL, INTENT(IN) :: zstd(klon) 
    REAL, INTENT(IN) :: zsig(klon) 
    REAL, INTENT(IN) :: zgam(klon) 
    REAL, INTENT(IN) :: zthe(klon)

    ALLOCATE(zmea_gw(klon), zpic_gw(klon), zval_gw(klon), zstd_gw(klon), zsig_gw(klon), zgam_gw(klon), zthe_gw(klon)) 
    
    zmea_gw(:)=zmea(:)
    zpic_gw(:)=zpic(:)
    zval_gw(:)=zval(:)
    zstd_gw(:)=zstd(:)
    zsig_gw(:)=zsig(:)
    zgam_gw(:)=zgam(:)
    zthe_gw(:)=zthe(:)

  END SUBROUTINE init_param_gw




  SUBROUTINE create_etat0_unstruct
  USE dimphy
  USE lmdz_xios
  USE infotrac_phy
  USE fonte_neige_mod
  USE pbl_surface_mod
  USE phys_state_var_mod
  USE indice_sol_mod
  USE surface_data,      ONLY: landice_opt
  USE mod_phys_lmdz_para
  USE print_control_mod, ONLY: lunout
  USE geometry_mod
  USE ioipsl_getin_p_mod, ONLY: getin_p

  IMPLICIT NONE
  INCLUDE 'dimsoil.h'
  include "clesphys.h"

    LOGICAL :: no_ter_antartique   ! If true, no land points are allowed at Antartic
    REAL,    DIMENSION(klon)                 :: tsol
    REAL,    DIMENSION(klon)                 :: sn
    REAL,    DIMENSION(klon)                 :: rugmer
    REAL,    DIMENSION(klon)                 :: run_off_lic_0
    REAL,    DIMENSION(klon)                 :: lic
    REAL,    DIMENSION(klon)                 :: fder

    REAL,    DIMENSION(klon,nbsrf)           :: qsurf, snsrf
    REAL,    DIMENSION(klon,nsoilmx,nbsrf)   :: tsoil
    
    REAL,    DIMENSION(klon_mpi)             :: tsol_mpi, qsol_mpi, zmasq_mpi, lic_mpi
    REAL,    DIMENSION(klon_mpi)             :: zmea_mpi, zstd_mpi, zsig_mpi, zgam_mpi, zthe_mpi
    REAL,    DIMENSION(klon_mpi)             :: cell_area_mpi
    REAL,    DIMENSION(klon_mpi,nbsrf)       :: pctsrf_mpi

    INCLUDE "compbl.h"
    INCLUDE "alpale.h"
   
    INTEGER :: ji,j,i
 

!--- Initial atmospheric CO2 conc. from .def file
    co2_ppm0 = co2_ppm

    IF (is_omp_master) THEN
      CALL xios_recv_field("ts",tsol_mpi)
      CALL xios_recv_field("qs",qsol_mpi)
      CALL xios_recv_field("mask",zmasq_mpi)
      IF (landice_opt .LT. 2) CALL xios_recv_field("landice",lic_mpi)
    ENDIF
    CALL scatter_omp(tsol_mpi,tsol)
    CALL scatter_omp(qsol_mpi,qsol)
    CALL scatter_omp(zmasq_mpi,zmasq)
    IF (landice_opt .LT. 2) CALL scatter_omp(lic_mpi,lic)

    radsol(:)   = 0.0
    rugmer(:) = 0.001
    sn(:)     = 0

    WHERE(qsol(:)<0) qsol(:)=0 
        
    WHERE(   zmasq(:)<EPSFRA) zmasq(:)=0.
    WHERE(1.-zmasq(:)<EPSFRA) zmasq(:)=1.

    pctsrf(:,:) = 0
    IF (landice_opt .LT. 2) THEN
       pctsrf(:,is_lic)=lic
       WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0.
       WHERE(zmasq(:)<EPSFRA)         pctsrf(:,is_lic)=0.

       pctsrf(:,is_ter)=zmasq(:)
       
       !--- Adequation with soil/sea mask
       DO ji=1,klon
          IF(zmasq(ji)>EPSFRA) THEN 
             IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN
                pctsrf(ji,is_lic)=zmasq(ji)
                pctsrf(ji,is_ter)=0.
             ELSE
                pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic)
                IF(pctsrf(ji,is_ter)<EPSFRA) THEN
                   pctsrf(ji,is_ter)=0.
                   pctsrf(ji,is_lic)=zmasq(ji)
                END IF
             END IF
          END IF
       END DO
   
    ELSE
       ! landice_opt=>2 : no land ice
       pctsrf(:,is_lic)=0.0
       pctsrf(:,is_ter)=zmasq(:)
    END IF





  !--- Option no_ter_antartique removes all land fractions souther than 60S.
  !--- Land ice is set instead of the land fractions on these latitudes.
  !--- The ocean and sea-ice fractions are not changed.
  !--- This option is only available if landice_opt<2.   
  IF (landice_opt .LT. 2) THEN
     no_ter_antartique=.FALSE.
     CALL getin_p('no_ter_antartique',no_ter_antartique)
     WRITE(lunout,*)"no_ter_antartique=",no_ter_antartique
     IF (no_ter_antartique) THEN
        ! Remove all land fractions souther than 60S and set land-ice instead
        WRITE(lunout,*) "Remove land fractions souther than 60deg south by increasing"
        WRITE(lunout,*) "the continental ice fractions. No land can now be found at Antartic."
        DO ji=1, klon
           IF (latitude_deg(ji)<-60.0) THEN
              pctsrf(ji,is_lic) = pctsrf(ji,is_lic) + pctsrf(ji,is_ter)
              pctsrf(ji,is_ter) = 0
           END IF
        END DO
     END IF
  END IF
    
! sub-surface ocean and sea ice (sea ice set to zero for start)
!*******************************************************************************
    pctsrf(:,is_oce)=(1.-zmasq(:))
    WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0.
    
!    zval(:)=max(0.,zmea-2*zstd(:))
!    zpic(:)=zmea+2*zstd(:)
    
!! WARNING    DON'T FORGET FOR LATER
!!ym  IF(couple) pctsrf(:,is_oce)=ocemask_fi(:)
!! 
    
! Init: tsol, qsol, sn, evap, tsoil, rain_fall, snow_fall, solsw, sollw, frugs
!*******************************************************************************
    DO i=1,nbsrf
     ftsol(:,i) = tsol
    END DO
  
    DO i=1,nbsrf
     snsrf(:,i) = sn
    END DO
!albedo SB >>>
!ym error : the sub surface dimension is the third not second
!    falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
!    falb_dir(:,is_oce,:)=0.5;  falb_dir(:,is_sic,:)=0.6
    falb_dir(:,:,is_ter)=0.08; falb_dir(:,:,is_lic)=0.6
    falb_dir(:,:,is_oce)=0.5;  falb_dir(:,:,is_sic)=0.6

!ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
!ym probably the uninitialized value was 0 for standard (regular grid) case
    falb_dif(:,:,:)=0
    u10m(:,:)=0  
    v10m(:,:)=0  
    treedrg(:,:,:)=0
!albedo SB <<<
    fevap(:,:) = 0.
    qsurf = 0.
    
    DO i=1,nbsrf 
      DO j=1,nsoilmx
        tsoil(:,j,i) = tsol
      END DO
    END DO
 
    rain_fall = 0.; snow_fall = 0.
    solsw = 165.;   sollw = -53.
    solswfdiff = 1. 
!ym warning missing init for sollwdown => set to 0
  sollwdown  = 0.
    
    
    t_ancien = 273.15
    u_ancien=0
    v_ancien=0
    q_ancien = 0.
    ql_ancien = 0.
    qs_ancien = 0.
    prlw_ancien = 0.
    prsw_ancien = 0.
    prw_ancien = 0.
    agesno = 0.

    wake_delta_pbl_TKE(:,:,:)=0
    wake_dens(:)=0
    awake_dens = 0.
    cv_gen = 0.
    ale_bl = 0.
    ale_bl_trig =0.
    alp_bl=0.
    ale_wake=0.
    ale_bl_stat=0.
    
    z0m(:,:)=0 ! ym missing 5th subsurface initialization
    z0m(:,is_oce) = rugmer(:)
    z0m(:,is_ter) = 0.01 ! MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
    z0m(:,is_lic) = 0.001 !MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
    z0m(:,is_sic) = 0.001
    z0h(:,:)=z0m(:,:)

    fder = 0.0
    clwcon = 0.0
    rnebcon = 0.0
    ratqs = 0.0
    run_off_lic_0 = 0.0 
    rugoro = 0.0

! Before phyredem calling, surface modules and values to be saved in startphy.nc
! are initialized
!*******************************************************************************
    pbl_tke(:,:,:) = 1.e-8 
    zmax0(:) = 40.
    f0(:) = 1.e-5
    sig1(:,:) = 0.
    w01(:,:) = 0.
    wake_deltat(:,:) = 0.
    wake_deltaq(:,:) = 0.
    wake_s(:) = 0.
    wake_cstar(:) = 0.
    wake_fip(:) = 0.
    wake_pe = 0.
    fm_therm = 0.
    entr_therm = 0.
    detr_therm = 0.
    awake_s = 0.

    CALL fonte_neige_init(run_off_lic_0)
    CALL pbl_surface_init( fder, snsrf, qsurf, tsoil )

    IF (iflag_pbl>1 .AND. iflag_wake>=1  .AND. iflag_pbl_split >=1) then
     delta_tsurf = 0.
     beta_aridity = 0.
    END IF
    ratqs_inter_ = 0.002
    rneb_ancien = 0.

    CALL gather_omp(cell_area,cell_area_mpi)
    CALL gather_omp(pctsrf,pctsrf_mpi)
    IF (is_omp_master) THEN
      CALL xios_send_field("area_ce0l",cell_area_mpi)
      CALL xios_send_field("fract_oce_ce0l",pctsrf_mpi(:,is_oce))
      CALL xios_send_field("fract_sic_ce0l",pctsrf_mpi(:,is_sic))
    ENDIF
    
    zmea(:) = zmea_gw(:)
    zpic(:) = zpic_gw(:)
    zval(:) = zval_gw(:)
    zstd(:) = zstd_gw(:)
    zsig(:) = zsig_gw(:)
    zgam(:) = zgam_gw(:)
    zthe(:) = zthe_gw(:)
    DEALLOCATE(zmea_gw, zpic_gw, zval_gw, zstd_gw, zsig_gw, zgam_gw, zthe_gw)

    CALL phyredem( "startphy.nc" )

  END SUBROUTINE create_etat0_unstruct


END MODULE create_etat0_unstruct_mod