Index: /LMDZ6/branches/DYNAMICO-conv/libf/phy_common/xios_writefield.F90 =================================================================== --- /LMDZ6/branches/DYNAMICO-conv/libf/phy_common/xios_writefield.F90 (revision 3323) +++ /LMDZ6/branches/DYNAMICO-conv/libf/phy_common/xios_writefield.F90 (revision 3323) @@ -0,0 +1,294 @@ +MODULE xios_writefield_mod +#ifdef CPP_XIOS + PRIVATE + + INTEGER,PARAMETER :: size_name=255 + INTEGER, PARAMETER :: MaxWriteField = 100 + CHARACTER(len=size_name), SAVE :: FieldName(MaxWriteField) +!$OMP THREADPRIVATE(FieldName) + INTEGER, SAVE :: FieldIt (MaxWriteField) +!$OMP THREADPRIVATE(FieldIt) + INTEGER,SAVE :: NbField = 0 +!$OMP THREADPRIVATE(NbField) + + INTERFACE xios_writefield + MODULE PROCEDURE xios_writefield2d,xios_writefield3d + END INTERFACE + + LOGICAL,SAVE :: output_native_grid +!$OMP THREADPRIVATE( output_native_grid) + + INTEGER,SAVE :: ni_glo , nj_glo +!$OMP THREADPRIVATE( ni_glo, nj_glo) + + + PUBLIC xios_writefield + +CONTAINS + + FUNCTION NameId(name_in) + IMPLICIT NONE + CHARACTER(LEN=*),INTENT(IN) :: name_in + INTEGER :: nameId + CHARACTER(LEN=size_name) :: name + INTEGER :: n + + name=name_in + + DO n=1,NbField + IF (name==fieldName(n)) THEN + nameId=n + RETURN + ENDIF + ENDDO + + nameId=0 + RETURN + END FUNCTION NameId + + + + SUBROUTINE xios_writefield2d(field,name_in) + USE dimphy + USE mod_phys_lmdz_para + USE xios + USE print_control_mod, ONLY: lunout + USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, klon_glo, grid_type, unstructured, regular_lonlat + USE wxios, ONLY: wxios_domain_param_unstructured, wxios_domain_param, wxios_set_context + USE ioipsl + IMPLICIT NONE + REAL, DIMENSION(:), INTENT(IN) :: field + CHARACTER(LEN=*) :: name_in + CHARACTER(LEN=size_name) :: name + TYPE(xios_context) :: xios_ctx + TYPE(xios_domaingroup) :: domain_definition + TYPE(xios_domaingroup) :: my_domaingroup + TYPE(xios_domain) :: my_domain + TYPE(xios_domain) :: my_regular_domain + TYPE(xios_filegroup) :: file_definition + TYPE(xios_file) :: my_file + TYPE(xios_field) :: my_field + TYPE(xios_fieldgroup) :: field_definition + TYPE(xios_generate_rectilinear_domain) :: generate_domain + TYPE(xios_interpolate_domain) :: interpolate_domain + INTEGER :: id + REAL,DIMENSION(klon_mpi) :: field_mpi + REAL :: Field2d(nbp_lon,jj_nb) + + + IF (size(field,1) /= klon_omp) THEN + WRITE(lunout,*) 'xios_writefield :: '//FieldName//' is not on the model grid' + RETURN + ENDIF + + name = TRIM(ADJUSTL(name_in)) + id=nameId(name) + + IF (id/=0) THEN + IF (is_omp_master) THEN + CALL xios_get_handle('context_lmdz_'//TRIM(name), xios_ctx) !Récupération + CALL xios_set_current_context(xios_ctx) !Activation + ENDIF + ELSE + output_native_grid=.FALSE. + ni_glo=0 + nj_glo=0 + IF (is_master) CALL getin("xios_writefield_nlon",ni_glo) + IF (is_master) CALL getin("xios_writefield_nlat",nj_glo) + CALL bcast(ni_glo) + CALL bcast(nj_glo) + IF (ni_glo==0 .OR. nj_glo==0) output_native_grid=.TRUE. + + IF (is_omp_master) THEN + CALL xios_context_initialize('context_lmdz_'//TRIM(name), COMM_LMDZ_PHY) + CALL xios_get_handle('context_lmdz_'//TRIM(name), xios_ctx) !Récupération + CALL xios_set_current_context(xios_ctx) !Activation + + CALL xios_define_calendar("D360") + CALL xios_set_start_date(xios_date(2000,1,1,0,0,0)) + CALL xios_set_time_origin(xios_date(2000,1,1,0,0,0)) + CALL xios_set_timestep(xios_second) + + CALL xios_get_handle("domain_definition",domain_definition) + CALL xios_add_child(domain_definition,my_domaingroup,"domaingroup") + CALL xios_add_child(my_domaingroup,my_domain,"domain") + IF (grid_type==unstructured .AND. .NOT. output_native_grid ) THEN + CALL xios_add_child(domain_definition, my_regular_domain, "regular_domain") + CALL xios_set_attr(my_regular_domain,ni_glo=ni_glo, nj_glo=nj_glo, type="rectilinear") + CALL xios_add_child(my_regular_domain,generate_domain) + CALL xios_set_attr(generate_domain,lon_start=-180., lat_start=90., lat_end=-90.) + CALL xios_add_child(my_regular_domain,interpolate_domain) + ENDIF + + ENDIF + + IF (grid_type==regular_lonlat) THEN + CALL wxios_domain_param("domain") + ELSE IF (grid_type==unstructured) THEN + CALL wxios_domain_param_unstructured("domaingroup") + ENDIF + + IF (is_omp_master) THEN + CALL xios_get_handle("file_definition",file_definition) + CALL xios_add_child(file_definition, my_file) + CALL xios_set_attr(my_file,name=TRIM(name),output_freq=xios_timestep,sync_freq=xios_timestep,type="one_file") + CALL xios_get_handle("field_definition",field_definition) + CALL xios_add_child(field_definition, my_field, TRIM(name)) + CALL xios_set_attr(my_field,domain_ref="domain",operation="instant") + CALL xios_add_child(my_file, my_field) + CALL xios_set_attr(my_field,field_ref=TRIM(name)) + IF (grid_type==unstructured .AND. .NOT. output_native_grid) CALL xios_set_attr(my_field,domain_ref="regular_domain") + + CALL xios_close_context_definition() + + NbField=NbField+1 + FieldName(NbField)=TRIM(name) + FieldIt(NbField)=0 + id=NbField + ENDIF + ENDIF + + CALL Gather_omp(field,field_mpi) + + IF (is_omp_master) THEN + FieldIt(id)=FieldIt(id)+1 + CALL xios_update_calendar(FieldIt(id)) + + IF (grid_type==regular_lonlat) THEN + CALL grid1Dto2D_mpi(field_mpi,Field2d) + CALL xios_send_field(TRIM(name), Field2d) + ELSE IF (grid_type==unstructured) THEN + CALL xios_send_field(TRIM(name), field_mpi) + ENDIF + ENDIF + + CALL wxios_set_context() + + END SUBROUTINE xios_writefield2d + + SUBROUTINE xios_writefield3d(field,name_in) + USE dimphy + USE mod_phys_lmdz_para + USE xios + USE print_control_mod, ONLY: lunout + USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, klon_glo, grid_type, unstructured, regular_lonlat + USE wxios, ONLY: wxios_domain_param_unstructured, wxios_domain_param, wxios_set_context + USE ioipsl + IMPLICIT NONE + REAL, DIMENSION(:,:), INTENT(IN) :: field + CHARACTER(LEN=*) :: name_in + CHARACTER(LEN=size_name) :: name + TYPE(xios_context) :: xios_ctx + TYPE(xios_domaingroup) :: domain_definition + TYPE(xios_domaingroup) :: my_domaingroup + TYPE(xios_domain) :: my_domain + TYPE(xios_domain) :: my_regular_domain + TYPE(xios_filegroup) :: file_definition + TYPE(xios_file) :: my_file + TYPE(xios_field) :: my_field + TYPE(xios_fieldgroup) :: field_definition + TYPE(xios_generate_rectilinear_domain) :: generate_domain + TYPE(xios_interpolate_domain) :: interpolate_domain + TYPE(xios_axisgroup) :: axis_definition + TYPE(xios_axis) :: my_axis + INTEGER :: id + REAL,DIMENSION(klon_mpi,size(field,2)) :: field_mpi + REAL :: Field2d(nbp_lon,jj_nb,size(field,2)) + REAL :: axis_value(size(field,2)) + INTEGER :: i + + IF (size(field,1) /= klon_omp) THEN + WRITE(lunout,*) 'xios_writefield :: '//FieldName//' is not on the model grid' + RETURN + ENDIF + + name = TRIM(ADJUSTL(name_in)) + id=nameId(name) + + IF (id/=0) THEN + + IF (is_omp_master) CALL xios_get_handle('context_lmdz_'//TRIM(name), xios_ctx) !Récupération + IF (is_omp_master) CALL xios_set_current_context(xios_ctx) !Activation + + ELSE + output_native_grid=.FALSE. + ni_glo=0 + nj_glo=0 + IF (is_master) CALL getin("xios_writefield_nlon",ni_glo) + IF (is_master) CALL getin("xios_writefield_nlat",nj_glo) + CALL bcast(ni_glo) + CALL bcast(nj_glo) + IF (ni_glo==0 .OR. nj_glo==0) output_native_grid=.TRUE. + + IF (is_omp_master) THEN + CALL xios_context_initialize('context_lmdz_'//TRIM(name), COMM_LMDZ_PHY) + CALL xios_get_handle('context_lmdz_'//TRIM(name), xios_ctx) !Récupération + CALL xios_set_current_context(xios_ctx) !Activation + + CALL xios_define_calendar("D360") + CALL xios_set_start_date(xios_date(2000,1,1,0,0,0)) + CALL xios_set_time_origin(xios_date(2000,1,1,0,0,0)) + CALL xios_set_timestep(xios_second) + + CALL xios_get_handle("domain_definition",domain_definition) + CALL xios_add_child(domain_definition,my_domaingroup,"domaingroup") + CALL xios_add_child(my_domaingroup,my_domain,"domain") + IF (grid_type==unstructured .AND. .NOT. output_native_grid) THEN + CALL xios_add_child(domain_definition, my_regular_domain, "regular_domain") + CALL xios_set_attr(my_regular_domain,ni_glo=ni_glo, nj_glo=nj_glo, type="rectilinear") + CALL xios_add_child(my_regular_domain,generate_domain) + CALL xios_set_attr(generate_domain,lon_start=-180., lat_start=90., lat_end=-90.) + CALL xios_add_child(my_regular_domain,interpolate_domain) + ENDIF + CALL xios_get_handle("axis_definition",axis_definition) + CALL xios_add_child(axis_definition,my_axis,"axis") + axis_value=(/(i,i=1,size(field,2))/) + CALL xios_set_attr(my_axis,name="z",n_glo=size(field,2),value=axis_value, unit="level" ) + ENDIF + + IF (grid_type==regular_lonlat) THEN + CALL wxios_domain_param("domain") + ELSE IF (grid_type==unstructured) THEN + CALL wxios_domain_param_unstructured("domaingroup") + ENDIF + + IF (is_omp_master) THEN + CALL xios_get_handle("file_definition",file_definition) + CALL xios_add_child(file_definition, my_file) + CALL xios_set_attr(my_file,name=TRIM(name),output_freq=xios_timestep,sync_freq=xios_timestep,type="one_file") + CALL xios_get_handle("field_definition",field_definition) + CALL xios_add_child(field_definition, my_field, TRIM(name)) + CALL xios_set_attr(my_field,domain_ref="domain",axis_ref="axis", operation="instant") + CALL xios_add_child(my_file, my_field) + CALL xios_set_attr(my_field,field_ref=TRIM(name)) + IF (grid_type==unstructured .AND. .NOT. output_native_grid) CALL xios_set_attr(my_field,domain_ref="regular_domain") + + CALL xios_close_context_definition() + + NbField=NbField+1 + FieldName(NbField)=TRIM(name) + FieldIt(NbField)=0 + id=NbField + ENDIF + ENDIF + + CALL Gather_omp(field,field_mpi) + + IF (is_omp_master) THEN + FieldIt(id)=FieldIt(id)+1 + CALL xios_update_calendar(FieldIt(id)) + + IF (grid_type==regular_lonlat) THEN + CALL grid1Dto2D_mpi(field_mpi,Field2d) + CALL xios_send_field(TRIM(name), Field2d) + ELSE IF (grid_type==unstructured) THEN + CALL xios_send_field(TRIM(name), field_mpi) + ENDIF + ENDIF + + CALL wxios_set_context() + + END SUBROUTINE xios_writefield3d + +#endif +END MODULE Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/create_etat0_limit_unstruct.F90 =================================================================== --- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/create_etat0_limit_unstruct.F90 (revision 3323) +++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/create_etat0_limit_unstruct.F90 (revision 3323) @@ -0,0 +1,98 @@ +MODULE etat0_limit_unstruct_mod + + LOGICAL, SAVE :: create_etat0_limit +!$OMP THREADPRIVATE(create_etat0_limit) + + + + +CONTAINS + + SUBROUTINE init_etat0_limit_unstruct + USE xios + USE mod_phys_lmdz_para + USE mod_grid_phy_lmdz + USE ioipsl, ONLY : getin, ioget_year_len + USE time_phylmdz_mod, ONLY : annee_ref + IMPLICIT NONE + + INTEGER :: iflag_phys,i + INTEGER :: ndays + REAL,ALLOCATABLE :: value(:) + + IF (grid_type==unstructured) THEN + IF (is_master) CALL getin("iflag_phys",iflag_phys) + CALL bcast(iflag_phys) + + IF (is_master) CALL getin('create_etat0_limit',create_etat0_limit) + CALL bcast(create_etat0_limit) + + ndays=ioget_year_len(annee_ref) + ALLOCATE(value(ndays)) + DO i=1,ndays + value(i)=i-1 + ENDDO + + IF (is_omp_master) CALL xios_set_axis_attr("time_year",n_glo=ndays,value=value) + + IF (create_etat0_limit) THEN + IF (iflag_phys<100) THEN + IF (is_omp_master) CALL xios_set_fieldgroup_attr("etat0_limit_read",read_access=.TRUE.,enabled=.TRUE.) + IF (is_omp_master) CALL xios_set_filegroup_attr("etat0_limit_read",enabled=.TRUE.) + ENDIF + IF (is_omp_master) CALL xios_set_file_attr("limit_write",enabled=.TRUE.) + ENDIF + + ENDIF + + END SUBROUTINE init_etat0_limit_unstruct + + SUBROUTINE create_etat0_limit_unstruct + USE mod_grid_phy_lmdz + USE create_etat0_unstruct_mod + USE create_limit_unstruct_mod + USE phyaqua_mod + USE mod_phys_lmdz_para + USE ioipsl, ONLY : getin + USE dimphy + USE xios + IMPLICIT NONE + INTEGER :: iflag_phys + + IF (grid_type==unstructured) THEN + + IF (is_master) CALL getin("iflag_phys",iflag_phys) + CALL bcast(iflag_phys) + + IF (iflag_phys<100) THEN + IF ( create_etat0_limit) THEN + CALL create_etat0_unstruct + CALL create_limit_unstruct + IF (is_omp_master) THEN + CALL xios_context_finalize() + CALL xios_set_current_context("icosagcm") ! very bad, need to find an other solution + CALL xios_context_finalize() + CALL xios_finalize() + ENDIF +!$OMP BARRIER + STOP 'create_etat0_limit_unstruct, Initial state file are created, all is fine' + ENDIF + ELSE + IF (create_etat0_limit) THEN + CALL iniaqua(klon,iflag_phys) + IF (is_omp_master) THEN + CALL xios_context_finalize() + CALL xios_set_current_context("icosagcm") ! very bad, need to find an other solution + CALL xios_context_finalize() + CALL xios_finalize() + ENDIF +!$OMP BARRIER + STOP 'create_etat0_limit_unstruct, Initial state file are created, all is fine' + ENDIF + ENDIF + ENDIF + + END SUBROUTINE create_etat0_limit_unstruct + +END MODULE etat0_limit_unstruct_mod + Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/create_etat0_unstruct.F90 =================================================================== --- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/create_etat0_unstruct.F90 (revision 3323) +++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/create_etat0_unstruct.F90 (revision 3323) @@ -0,0 +1,186 @@ +MODULE create_etat0_unstruct_mod + + + + + + +CONTAINS + + SUBROUTINE create_etat0_unstruct + USE dimphy +#ifdef CPP_XIOS + USE xios + USE infotrac_phy + USE fonte_neige_mod + USE pbl_surface_mod + USE phys_state_var_mod + USE indice_sol_mod + USE mod_phys_lmdz_para + IMPLICIT NONE + INCLUDE 'dimsoil.h' + + 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) :: qsolsrf, 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 + + INTEGER :: ji,j,i + + 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) + CALL xios_recv_field("landice",lic_mpi) + CALL xios_recv_field("zmea",zmea_mpi) + CALL xios_recv_field("zstd",zstd_mpi) + CALL xios_recv_field("zsig",zsig_mpi) + CALL xios_recv_field("zgam",zgam_mpi) + CALL xios_recv_field("zthe",zthe_mpi) + ENDIF + CALL scatter_omp(tsol_mpi,tsol) + CALL scatter_omp(qsol_mpi,qsol) + CALL scatter_omp(zmasq_mpi,zmasq) + CALL scatter_omp(lic_mpi,lic) + CALL scatter_omp(zmea_mpi,zmea) + CALL scatter_omp(zstd_mpi,zstd) + CALL scatter_omp(zsig_mpi,zsig) + CALL scatter_omp(zgam_mpi,zgam) + CALL scatter_omp(zthe_mpi,zthe) + + radsol(:) = 0.0 + rugmer(:) = 0.001 + sn(:) = 0 + + WHERE(qsol(:)<0) qsol(:)=0 + + WHERE( zmasq(:)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)>> +!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 + +!albedo SB <<< + fevap(:,:) = 0. + DO i=1,nbsrf + qsolsrf(:,i)=150. + END DO + + 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. + t_ancien = 273.15 + u_ancien=0 + v_ancien=0 + q_ancien = 0. + agesno = 0. + + z0m(:,is_oce) = rugmer(:) + + z0m(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0) + z0m(:,is_lic) = 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. + ale_bl = 0. + ale_bl_trig =0. + alp_bl =0. + CALL fonte_neige_init(run_off_lic_0) + CALL pbl_surface_init( fder, snsrf, qsolsrf, tsoil ) + CALL phyredem( "startphy.nc" ) + +#endif + END SUBROUTINE create_etat0_unstruct + + +END MODULE create_etat0_unstruct_mod Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/create_limit_unstruct.F90 =================================================================== --- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/create_limit_unstruct.F90 (revision 3323) +++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/create_limit_unstruct.F90 (revision 3323) @@ -0,0 +1,277 @@ +MODULE create_limit_unstruct_mod + PRIVATE + INTEGER, PARAMETER :: lmdep=12 + + PUBLIC create_limit_unstruct + +CONTAINS + + SUBROUTINE create_limit_unstruct + USE dimphy +#ifdef CPP_XIOS + USE xios + USE ioipsl, ONLY : ioget_year_len + USE time_phylmdz_mod, ONLY : annee_ref + USE indice_sol_mod + USE phys_state_var_mod + USE mod_phys_lmdz_para + IMPLICIT NONE + INCLUDE "iniprint.h" + REAL, DIMENSION(klon,lmdep) :: sic + REAL, DIMENSION(klon,lmdep) :: sst + REAL, DIMENSION(klon,lmdep) :: rugos + REAL, DIMENSION(klon,lmdep) :: albedo + REAL, DIMENSION(klon_mpi,lmdep) :: sic_mpi + REAL, DIMENSION(klon_mpi,lmdep) :: sst_mpi + REAL, DIMENSION(klon_mpi,lmdep) :: rugos_mpi + REAL, DIMENSION(klon_mpi,lmdep) :: albedo_mpi + INTEGER :: ndays + REAL :: fi_ice(klon) + REAL, ALLOCATABLE :: sic_year(:,:) + REAL, ALLOCATABLE :: sst_year(:,:) + REAL, ALLOCATABLE :: rugos_year(:,:) + REAL, ALLOCATABLE :: albedo_year(:,:) + REAL, ALLOCATABLE :: pctsrf_t(:,:,:) + REAL, ALLOCATABLE :: phy_bil(:,:) + REAL, ALLOCATABLE :: sst_year_mpi(:,:) + REAL, ALLOCATABLE :: rugos_year_mpi(:,:) + REAL, ALLOCATABLE :: albedo_year_mpi(:,:) + REAL, ALLOCATABLE :: pctsrf_t_mpi(:,:,:) + REAL, ALLOCATABLE :: phy_bil_mpi(:,:) + INTEGER :: l,k + INTEGER :: nbad + + ndays=ioget_year_len(annee_ref) + + IF (is_omp_master) THEN + CALL xios_recv_field("sic_limit",sic_mpi) + CALL xios_recv_field("sst_limit",sst_mpi) + CALL xios_recv_field("rugos_limit",rugos_mpi) + CALL xios_recv_field("albedo_limit",albedo_mpi) + ENDIF + CALL scatter_omp(sic_mpi,sic) + CALL scatter_omp(sst_mpi,sst) + CALL scatter_omp(rugos_mpi,rugos) + CALL scatter_omp(albedo_mpi,albedo) + + ALLOCATE(sic_year(klon,ndays)) + ALLOCATE(sst_year(klon,ndays)) + ALLOCATE(rugos_year(klon,ndays)) + ALLOCATE(albedo_year(klon,ndays)) + ALLOCATE(pctsrf_t(klon,nbsrf,ndays)) + ALLOCATE(phy_bil(klon,ndays)); phy_bil=0.0 + +! sic + CALL time_interpolation(ndays,sic,'gregorian',sic_year) + sic_year(:,:)=sic_year(:,:)/100. ! convert percent to fraction + WHERE(sic_year(:,:)>1.0) sic_year(:,:)=1.0 ! Some fractions have some time large negative values + WHERE(sic_year(:,:)<0.0) sic_year(:,:)=0.0 ! probably better to apply alse this filter before horizontal interpolation + +! sst + CALL time_interpolation(ndays,sst,'gregorian',sst_year) + WHERE(sst_year(:,:)<271.38) sst_year(:,:)=271.38 + +! rugos + DO l=1, lmdep + WHERE(NINT(zmasq(:))/=1) rugos(:,l)=0.001 + ENDDO + CALL time_interpolation(ndays,rugos,'360d',rugos_year) + +! albedo + CALL time_interpolation(ndays,albedo,'360d',albedo_year) + + + DO k=1,ndays + fi_ice=sic_year(:,k) + WHERE(fi_ice>=1.0 ) fi_ice=1.0 + WHERE(fi_ice=1.0-zmasq) + pctsrf_t(:,is_sic,k)=1.0-zmasq + pctsrf_t(:,is_oce,k)=0.0 + ELSEWHERE + pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k) + WHERE(pctsrf_t(:,is_oce,k)0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad + nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA) + IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad + END DO + + ALLOCATE(sst_year_mpi(klon_mpi,ndays)) + ALLOCATE(rugos_year_mpi(klon_mpi,ndays)) + ALLOCATE(albedo_year_mpi(klon_mpi,ndays)) + ALLOCATE(pctsrf_t_mpi(klon_mpi,nbsrf,ndays)) + ALLOCATE(phy_bil_mpi(klon_mpi,ndays)) + + CALL gather_omp(pctsrf_t , pctsrf_t_mpi) + CALL gather_omp(sst_year , sst_year_mpi) + CALL gather_omp(phy_bil , phy_bil_mpi) + CALL gather_omp(albedo_year, albedo_year_mpi) + CALL gather_omp(rugos_year , rugos_year_mpi) + + IF (is_omp_master) THEN + CALL xios_send_field("foce_limout",pctsrf_t_mpi(:,is_oce,:)) + CALL xios_send_field("fsic_limout",pctsrf_t_mpi(:,is_sic,:)) + CALL xios_send_field("fter_limout",pctsrf_t_mpi(:,is_ter,:)) + CALL xios_send_field("flic_limout",pctsrf_t_mpi(:,is_lic,:)) + CALL xios_send_field("sst_limout", sst_year_mpi) + CALL xios_send_field("bils_limout",phy_bil_mpi) + CALL xios_send_field("alb_limout", albedo_year_mpi) + CALL xios_send_field("rug_limout", rugos_year_mpi) + ENDIF +#endif + END SUBROUTINE create_limit_unstruct + + + SUBROUTINE time_interpolation(ndays,field_in,calendar,field_out) + USE pchsp_95_m, only: pchsp_95 + USE pchfe_95_m, only: pchfe_95 + USE arth_m, only: arth + USE dimphy, ONLY : klon + USE ioipsl, ONLY : ioget_year_len + USE time_phylmdz_mod, ONLY : annee_ref + USE mod_phys_lmdz_para + IMPLICIT NONE + INCLUDE "iniprint.h" + + INTEGER, INTENT(IN) :: ndays + REAL, INTENT(IN) :: field_in(klon,lmdep) + CHARACTER(LEN=*),INTENT(IN) :: calendar + REAL, INTENT(OUT) :: field_out(klon,ndays) + + INTEGER :: ndays_in + REAL :: timeyear(lmdep) + REAL :: yder(lmdep) + INTEGER :: ij,ierr, n_extrap + LOGICAL :: skip + + + IF (is_omp_master) ndays_in=year_len(annee_ref, calendar) + CALL bcast_omp(ndays_in) + timeyear=mid_months(annee_ref, calendar, lmdep) + + n_extrap = 0 + skip=.FALSE. + DO ij=1,klon + yder = pchsp_95(timeyear, field_in(ij, :), ibeg=2, iend=2, vc_beg=0., vc_end=0.) + CALL pchfe_95(timeyear, field_in(ij, :), yder, skip, arth(0., real(ndays_in) / ndays, ndays), field_out(ij, :), ierr) + if (ierr < 0) stop 1 + n_extrap = n_extrap + ierr + END DO + + IF (n_extrap /= 0) then + WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap + ENDIF + + + END SUBROUTINE time_interpolation + !------------------------------------------------------------------------------- + ! + FUNCTION year_len(y,cal_in) + ! + !------------------------------------------------------------------------------- + USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len + IMPLICIT NONE + !------------------------------------------------------------------------------- + ! Arguments: + INTEGER :: year_len + INTEGER, INTENT(IN) :: y + CHARACTER(LEN=*), INTENT(IN) :: cal_in + !------------------------------------------------------------------------------- + ! Local variables: + CHARACTER(LEN=20) :: cal_out ! calendar (for outputs) + !------------------------------------------------------------------------------- + !--- Getting the input calendar to reset at the end of the function + CALL ioget_calendar(cal_out) + + !--- Unlocking calendar and setting it to wanted one + CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in)) + + !--- Getting the number of days in this year + year_len=ioget_year_len(y) + + !--- Back to original calendar + CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out)) + + END FUNCTION year_len + ! + !------------------------------------------------------------------------------- + + + !------------------------------------------------------------------------------- + ! + FUNCTION mid_months(y,cal_in,nm) + ! + !------------------------------------------------------------------------------- + USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len + IMPLICIT NONE + !------------------------------------------------------------------------------- + ! Arguments: + INTEGER, INTENT(IN) :: y ! year + CHARACTER(LEN=*), INTENT(IN) :: cal_in ! calendar + INTEGER, INTENT(IN) :: nm ! months/year number + REAL, DIMENSION(nm) :: mid_months ! mid-month times + !------------------------------------------------------------------------------- + ! Local variables: + CHARACTER(LEN=99) :: mess ! error message + CHARACTER(LEN=20) :: cal_out ! calendar (for outputs) + INTEGER, DIMENSION(nm) :: mnth ! months lengths (days) + INTEGER :: m ! months counter + INTEGER :: nd ! number of days + INTEGER :: k + !------------------------------------------------------------------------------- + nd=year_len(y,cal_in) + + IF(nm==12) THEN + + !--- Getting the input calendar to reset at the end of the function + CALL ioget_calendar(cal_out) + + !--- Unlocking calendar and setting it to wanted one + CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in)) + + !--- Getting the length of each month + DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO + + !--- Back to original calendar + CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out)) + + ELSE IF(MODULO(nd,nm)/=0) THEN + WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',& + nm,' months/year. Months number should divide days number.' + CALL abort_physic('mid_months',TRIM(mess),1) + + ELSE + mnth=(/(m,m=1,nm,nd/nm)/) + END IF + + !--- Mid-months times + mid_months(1)=0.5*REAL(mnth(1)) + DO k=2,nm + mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k)) + END DO + + END FUNCTION mid_months + + +END MODULE create_limit_unstruct_mod