Index: LMDZ6/trunk/libf/phylmd/create_etat0_unstruct.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/create_etat0_unstruct.F90 (revision 4620) +++ (revision ) @@ -1,270 +1,0 @@ -MODULE create_etat0_unstruct_mod - - - - - - -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 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' - - 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) :: 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 - REAL, DIMENSION(klon_mpi) :: cell_area_mpi - REAL, DIMENSION(klon_mpi,nbsrf) :: pctsrf_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) - IF (landice_opt .LT. 2) 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) - IF (landice_opt .LT. 2) 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)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)>> -!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. -!ym warning missing init for sollwdown => set to 0 - sollwdown = 0. - - - t_ancien = 273.15 - u_ancien=0 - v_ancien=0 - q_ancien = 0. - agesno = 0. - - 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. - 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 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 - - CALL phyredem( "startphy.nc" ) - - END SUBROUTINE create_etat0_unstruct - - -END MODULE create_etat0_unstruct_mod Index: LMDZ6/trunk/libf/phylmd/create_etat0_unstruct_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/create_etat0_unstruct_mod.F90 (revision 4621) +++ LMDZ6/trunk/libf/phylmd/create_etat0_unstruct_mod.F90 (revision 4621) @@ -0,0 +1,270 @@ +MODULE create_etat0_unstruct_mod + + + + + + +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 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' + + 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) :: 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 + REAL, DIMENSION(klon_mpi) :: cell_area_mpi + REAL, DIMENSION(klon_mpi,nbsrf) :: pctsrf_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) + IF (landice_opt .LT. 2) 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) + IF (landice_opt .LT. 2) 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)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)>> +!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. +!ym warning missing init for sollwdown => set to 0 + sollwdown = 0. + + + t_ancien = 273.15 + u_ancien=0 + v_ancien=0 + q_ancien = 0. + agesno = 0. + + 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. + 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 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 + + CALL phyredem( "startphy.nc" ) + + END SUBROUTINE create_etat0_unstruct + + +END MODULE create_etat0_unstruct_mod Index: LMDZ6/trunk/libf/phylmd/create_limit_unstruct.F90create_limit_unstruct_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/create_limit_unstruct.F90create_limit_unstruct_mod.F90 (revision 4620) +++ (revision ) @@ -1,317 +1,0 @@ -MODULE create_limit_unstruct_mod - PRIVATE - INTEGER, PARAMETER :: lmdep=12 - - PUBLIC create_limit_unstruct - -CONTAINS - - - SUBROUTINE create_limit_unstruct - USE dimphy - USE lmdz_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(:,:),ALLOCATABLE :: sic - REAL, DIMENSION(:,:),ALLOCATABLE :: sst - REAL, DIMENSION(klon,lmdep) :: rugos - REAL, DIMENSION(klon,lmdep) :: albedo - REAL, DIMENSION(:,:),ALLOCATABLE :: sic_mpi - REAL, DIMENSION(:,:),ALLOCATABLE :: 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 - INTEGER :: sic_time_axis_size - INTEGER :: sst_time_axis_size - CHARACTER(LEN=99) :: mess ! error message - - - ndays=ioget_year_len(annee_ref) - - IF (is_omp_master) CALL xios_get_axis_attr("time_sic",n_glo=sic_time_axis_size) - CALL bcast_omp(sic_time_axis_size) - ALLOCATE(sic_mpi(klon_mpi,sic_time_axis_size)) - ALLOCATE(sic(klon,sic_time_axis_size)) - - - IF (is_omp_master) CALL xios_get_axis_attr("time_sst",n_glo=sst_time_axis_size) - CALL bcast_omp(sst_time_axis_size) - ALLOCATE(sst_mpi(klon_mpi,sst_time_axis_size)) - ALLOCATE(sst(klon,sst_time_axis_size)) - - 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 - IF (sic_time_axis_size==lmdep) THEN - CALL time_interpolation(ndays,sic,'gregorian',sic_year) - ELSE IF (sic_time_axis_size==ndays) THEN - sic_year=sic - ELSE - WRITE(mess,*) 'sic time axis is nor montly, nor daily. sic time interpolation ',& - 'is requiered but is not currently managed' - CALL abort_physic('create_limit_unstruct',TRIM(mess),1) - ENDIF - - 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 - IF (sst_time_axis_size==lmdep) THEN - CALL time_interpolation(ndays,sst,'gregorian',sst_year) - ELSE IF (sst_time_axis_size==ndays) THEN - sst_year=sst - ELSE - WRITE(mess,*)'sic time axis is nor montly, nor daily. sic time interpolation ',& - 'is requiered but is not currently managed' - CALL abort_physic('create_limit_unstruct',TRIM(mess),1) - ENDIF - 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,'360_day',rugos_year) - -! albedo - CALL time_interpolation(ndays,albedo,'360_day',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 - 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 - - CHARACTER (len = 50) :: modname = 'create_limit_unstruct.time_interpolation' - CHARACTER (len = 80) :: abort_message - - - IF (is_omp_master) ndays_in=year_len(annee_ref, calendar) - CALL bcast_omp(ndays_in) - IF (is_omp_master) timeyear=mid_months(annee_ref, calendar, lmdep) - CALL bcast_omp(timeyear) - - 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) then - abort_message='error in pchfe_95' - CALL abort_physic(modname,abort_message,1) - endif - 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 Index: LMDZ6/trunk/libf/phylmd/create_limit_unstruct_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/create_limit_unstruct_mod.F90 (revision 4621) +++ LMDZ6/trunk/libf/phylmd/create_limit_unstruct_mod.F90 (revision 4621) @@ -0,0 +1,317 @@ +MODULE create_limit_unstruct_mod + PRIVATE + INTEGER, PARAMETER :: lmdep=12 + + PUBLIC create_limit_unstruct + +CONTAINS + + + SUBROUTINE create_limit_unstruct + USE dimphy + USE lmdz_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(:,:),ALLOCATABLE :: sic + REAL, DIMENSION(:,:),ALLOCATABLE :: sst + REAL, DIMENSION(klon,lmdep) :: rugos + REAL, DIMENSION(klon,lmdep) :: albedo + REAL, DIMENSION(:,:),ALLOCATABLE :: sic_mpi + REAL, DIMENSION(:,:),ALLOCATABLE :: 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 + INTEGER :: sic_time_axis_size + INTEGER :: sst_time_axis_size + CHARACTER(LEN=99) :: mess ! error message + + + ndays=ioget_year_len(annee_ref) + + IF (is_omp_master) CALL xios_get_axis_attr("time_sic",n_glo=sic_time_axis_size) + CALL bcast_omp(sic_time_axis_size) + ALLOCATE(sic_mpi(klon_mpi,sic_time_axis_size)) + ALLOCATE(sic(klon,sic_time_axis_size)) + + + IF (is_omp_master) CALL xios_get_axis_attr("time_sst",n_glo=sst_time_axis_size) + CALL bcast_omp(sst_time_axis_size) + ALLOCATE(sst_mpi(klon_mpi,sst_time_axis_size)) + ALLOCATE(sst(klon,sst_time_axis_size)) + + 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 + IF (sic_time_axis_size==lmdep) THEN + CALL time_interpolation(ndays,sic,'gregorian',sic_year) + ELSE IF (sic_time_axis_size==ndays) THEN + sic_year=sic + ELSE + WRITE(mess,*) 'sic time axis is nor montly, nor daily. sic time interpolation ',& + 'is requiered but is not currently managed' + CALL abort_physic('create_limit_unstruct',TRIM(mess),1) + ENDIF + + 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 + IF (sst_time_axis_size==lmdep) THEN + CALL time_interpolation(ndays,sst,'gregorian',sst_year) + ELSE IF (sst_time_axis_size==ndays) THEN + sst_year=sst + ELSE + WRITE(mess,*)'sic time axis is nor montly, nor daily. sic time interpolation ',& + 'is requiered but is not currently managed' + CALL abort_physic('create_limit_unstruct',TRIM(mess),1) + ENDIF + 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,'360_day',rugos_year) + +! albedo + CALL time_interpolation(ndays,albedo,'360_day',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 + 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 + + CHARACTER (len = 50) :: modname = 'create_limit_unstruct.time_interpolation' + CHARACTER (len = 80) :: abort_message + + + IF (is_omp_master) ndays_in=year_len(annee_ref, calendar) + CALL bcast_omp(ndays_in) + IF (is_omp_master) timeyear=mid_months(annee_ref, calendar, lmdep) + CALL bcast_omp(timeyear) + + 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) then + abort_message='error in pchfe_95' + CALL abort_physic(modname,abort_message,1) + endif + 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