Ignore:
Timestamp:
2024-01-14T15:15:57+01:00 (6 months ago)
Author:
jan.polcher
Message:

The modifications performed on the trunk for the coupling to WRF and the interpolation by XIOS on curviliean grids has been backported.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/ORCHIDEE_2_2/ORCHIDEE/src_oasisdriver/orchideeoasis.f90

    r7265 r8377  
    2222!! 
    2323!! SVN          : 
    24 !! $HeadURL: $ 
    25 !! $Date: $ 
    26 !! $Revision: $ 
     24!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/trunk/ORCHIDEE/src_oasisdriver/orchideeoasis.f90 $ 
     25!! $Date: 2023-08-17 23:22:02 +0200 (Thu, 17 Aug 2023) $ 
     26!! $Revision: 8140 $ 
    2727!! \n 
    2828!_ ================================================================================================================================ 
    … …  
    3535  USE netcdf 
    3636  ! 
     37  USE xios 
    3738  ! 
    3839  USE ioipsl_para 
    … …  
    4041  ! 
    4142  USE grid 
     43  USE time 
    4244  USE timer 
    43  
     45  USE constantes 
     46  USE constantes_soil 
    4447  USE globgrd 
    4548  USE orchoasis_tools 
    … …  
    4851  USE control 
    4952  USE ioipslctrl 
     53  USE xios_orchidee 
    5054  ! 
    5155  USE mod_oasis 
    … …  
    6569  REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: corners_glo 
    6670  REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: corners_lon, corners_lat 
     71  REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: corners_lon_tmp, corners_lat_tmp 
    6772  INTEGER(i_std) :: nbindex_g, kjpindex 
    6873  INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: kindex, kindex_g 
    … …  
    7883  !- 
    7984  INTEGER(i_std) :: i, j, ik, nbdt, first_point 
     85  INTEGER :: ierr 
    8086  INTEGER(i_std) :: itau, itau_offset, itau_sechiba 
    81   REAL(r_std)    :: date0, date0_shifted, dt, julian, julian0 
     87  REAL(r_std)    :: date0, date0_shifted, dt, julian 
    8288  INTEGER(i_std) :: rest_id, rest_id_stom 
    8389  INTEGER(i_std) ::  hist_id, hist2_id, hist_id_stom, hist_id_stom_IPCC 
    8490  REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lalo_loc 
    85   INTEGER(i_std) :: iim, jjm 
     91  INTEGER(i_std) :: iim, jjm, in, jk 
    8692  REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lon, lat 
    8793  !- 
    … …  
    110116  !- output fields 
    111117  !- 
    112   REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: z0            !! Surface roughness 
     118  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: z0m           !! Surface momentum roughness 
     119  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: z0h           !! Surface scalar transport roughness 
    113120  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: coastalflow   !! Diffuse flow of water into the ocean (m^3/dt) 
    114121  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: riverflow     !! Largest rivers flowing into the ocean (m^3/dt) 
    … …  
    128135  REAL(r_std), ALLOCATABLE, DIMENSION (:,:)          :: lai           !! Leaf area index (m^2 m^{-2} 
    129136  REAL(r_std), ALLOCATABLE, DIMENSION (:,:)          :: height        !! Vegetation Height (m) 
     137  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: run_off_lic      !! Calving, melting and liquid precipitation, needed when coupled to atmosphere 
     138  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: run_off_lic_frac !! Cell fraction corresponding to run_off_lic 
     139 
    130140  !- 
    131141  !- Declarations for OASIS 
    … …  
    178188  ! Set parallel processing in ORCHIDEE 
    179189  ! 
    180   CALL Init_orchidee_para(LOCAL_OASIS_COMM)  
     190  CALL Init_orchidee_para(LOCAL_OASIS_COMM, .TRUE.)  
    181191  !==================================================================================== 
    182192  ! 
    … …  
    275285  !- 
    276286  CALL orchoasis_time(date0, dt, nbdt) 
     287  CALL time_initialize(0, date0, dt, "END") 
    277288  !- 
    278289  !- lalo needs to be created before going into the parallel region 
    … …  
    355366  ! 
    356367  ALLOCATE(lon(iim,jjm), lat(iim,jjm)) 
     368  ALLOCATE(corners_lon(nbseg,iim,jj_para_nb(mpi_rank)), corners_lat(nbseg,iim,jj_para_nb(mpi_rank))) 
    357369  ALLOCATE(kindex(kjpindex)) 
    358370  ! 
    359371  lon=lon_glo(:,jj_para_begin(mpi_rank):jj_para_end(mpi_rank)) 
    360372  lat=lat_glo(:,jj_para_begin(mpi_rank):jj_para_end(mpi_rank)) 
     373  DO in=1,nbseg 
     374     DO j=1,jj_para_nb(mpi_rank) 
     375        jk=jj_para_begin(mpi_rank)+(j-1) 
     376        DO i=1,iim 
     377           ! Copy the corners/bounds of the meshes which are local to this processor. 
     378           ! Serves for OASIS and XIOS. 
     379           corners_lon(in,i,j)=corners_glo(i,jk,in,1) 
     380           corners_lat(in,i,j)=corners_glo(i,jk,in,2) 
     381        ENDDO 
     382     ENDDO 
     383  ENDDO 
    361384  ! 
    362385  ! 
    … …  
    407430  ALLOCATE(carbwh(kjpindex), carbha(kjpindex)) 
    408431  ALLOCATE(tsol_rad(kjpindex), temp_sol_new(kjpindex), qsurf(kjpindex)) 
    409   ALLOCATE(albedo(kjpindex,2), emis(kjpindex), z0(kjpindex)) 
     432  ALLOCATE(albedo(kjpindex,2), emis(kjpindex), z0m(kjpindex), z0h(kjpindex)) 
    410433  ALLOCATE(veget(kjpindex,nvm)) 
    411434  ALLOCATE(lai(kjpindex,nvm)) 
    412435  ALLOCATE(height(kjpindex,nvm)) 
     436  ALLOCATE(run_off_lic(kjpindex)) 
     437  ALLOCATE(run_off_lic_frac(kjpindex)) 
     438  
    413439  ! 
    414440  WRITE(numout,*) "Rank", mpi_rank, "Domain size per proc !:", iim_glo, jjm_glo, kjpindex, SIZE(kindex) 
    415441  WRITE(numout,*) "Rank", mpi_rank, "In parallel region land index starts at : ", kindex(1) 
    416   !  
    417   CALL orchoasis_time(date0, dt, nbdt) 
    418   ! 
    419   CALL control_initialize 
    420   dt_sechiba = dt 
    421   !  
    422   itau = 0 
    423   ! 
    424   CALL ioipslctrl_restini(itau, date0, dt, rest_id, rest_id_stom, itau_offset, date0_shifted) 
    425   ! 
    426   ! Get the date of the first time step 
    427   ! 
    428   julian = date0 + 0.5*(dt/one_day) 
    429   CALL ju2ymds (julian, year, month, day, sec) 
    430   ! 
    431   in_julian = itau2date(itau, date0, dt) 
    432   CALL ymds2ju (year,1,1,zero, julian0) 
    433   julian_diff = in_julian-julian0 
    434  
    435   CALL xios_orchidee_init( MPI_COMM_ORCH,                    & 
    436        date0,    year,    month,           day, julian_diff, & 
    437        lon,      lat ) 
    438   ! 
    439   itau_sechiba = itau + itau_offset 
    440   ! 
    441   WRITE(numout,*) "orchoasis_history :", iim, jjm, SHAPE(lon) 
    442   !- Initialize IOIPSL sechiba output files 
    443   CALL ioipslctrl_history(iim, jjm, lon, lat,  kindex, kjpindex, itau_sechiba, & 
    444        date0, dt, hist_id, hist2_id, hist_id_stom, hist_id_stom_IPCC) 
    445   WRITE(numout,*) "HISTORY : Define for ", itau, date0, dt 
    446442  ! 
    447443  !- 
    … …  
    462458     CALL oasis_write_grid(oasis_gridname, iim_glo, jjm_glo, lon_glo, lat_glo) 
    463459     ! 
    464      ALLOCATE(corners_lon(iim_glo, jjm_glo, nbseg), corners_lat(iim_glo, jjm_glo, nbseg)) 
    465      corners_lon(:,:,:) = corners_glo(:,:,:,1) 
    466      corners_lat(:,:,:) = corners_glo(:,:,:,2) 
    467      CALL oasis_write_corner(oasis_gridname, iim_glo, jjm_glo, nbseg, corners_lon, corners_lat) 
     460     ALLOCATE(corners_lon_tmp(iim_glo, jjm_glo, nbseg), corners_lat_tmp(iim_glo, jjm_glo, nbseg)) 
     461     corners_lon_tmp(:,:,:) = corners_glo(:,:,:,1) 
     462     corners_lat_tmp(:,:,:) = corners_glo(:,:,:,2) 
     463     CALL oasis_write_corner(oasis_gridname, iim_glo, jjm_glo, nbseg, corners_lon_tmp, corners_lat_tmp) 
    468464     ! 
    469465     ALLOCATE(maskinv_glo(iim_glo, jjm_glo)) 
    … …  
    490486        CALL oasis_write_grid(oasis_gridname, iim_glo, jjm_glo, lon_glo, lat_glo) 
    491487        ! 
    492         CALL oasis_write_corner(oasis_gridname, iim_glo, jjm_glo, nbseg, corners_lon, corners_lat) 
     488        CALL oasis_write_corner(oasis_gridname, iim_glo, jjm_glo, nbseg, corners_lon_tmp, corners_lat_tmp) 
    493489        ! 
    494490        CALL oasis_write_mask(oasis_gridname, iim_glo, jjm_glo, maskinv_glo) 
    … …  
    498494        CALL oasis_terminate_grids_writing() 
    499495     ENDIF 
    500      DEALLOCATE(corners_lon, corners_lat) 
     496     DEALLOCATE(corners_lon_tmp, corners_lat_tmp) 
    501497     ! 
    502498  ENDIF 
    … …  
    512508  CALL orchoasis_defvar(mpi_rank, kjpindex) 
    513509  ! 
     510  !--------------------------------------------------------------------------------------- 
     511  !- 
     512  !- Initialize the restart process and XIOS. 
     513  !- 
     514  !--------------------------------------------------------------------------------------- 
     515  CALL control_initialize 
     516  dt_sechiba = dt 
     517  !  
     518  itau = 0 
     519  ! 
     520  CALL ioipslctrl_restini(itau, date0, dt, rest_id, rest_id_stom, itau_offset, date0_shifted) 
     521  ! 
     522  ! Get the date of the first time step 
     523  ! 
     524  ! To ensure that itau starts with 0 at date0 for the restart, we have to set an off-set to achieve this.  
     525  ! itau_offset will get used to prduce itau_sechiba. 
     526  ! 
     527  itau_offset=-itau_offset-1 
     528  ! 
     529  ! Get the date of the first time step 
     530  ! 
     531  WRITE(numout,*) "itau_offset : itau_offset : ", itau_offset 
     532  WRITE(numout,*) "itau_offset : date0 ymd : ", year_end, month_end, day_end, sec_end 
     533  WRITE(numout,*) "itau_offset : date0, sec, dt : ", sec_end, dt 
     534   
     535  CALL xios_orchidee_init( MPI_COMM_ORCH,                    & 
     536       date0,    year_end, month_end,  day_end, julian_diff, & 
     537       lon,      lat,      corners_lon, corners_lat ) 
     538  ! 
     539  CALL sechiba_xios_initialize  
     540  ! 
     541  CALL xios_orchidee_close_definition 
     542  WRITE(numout,*) 'After xios_orchidee_close_definition' 
     543  ! 
     544  itau_sechiba = itau + itau_offset 
     545  ! 
     546  WRITE(numout,*) "orchoasis_history :", iim, jjm, SHAPE(lon) 
     547  !- Initialize IOIPSL sechiba output files 
     548  CALL ioipslctrl_history(iim, jjm, lon, lat,  kindex, kjpindex, itau_sechiba, & 
     549       date0, dt, hist_id, hist2_id, hist_id_stom, hist_id_stom_IPCC) 
     550  WRITE(numout,*) "HISTORY : Define for ", itau, date0, dt 
     551  ! 
    514552  !- 
    515553  !--------------------------------------------------------------------------------------- 
    … …  
    519557  !--------------------------------------------------------------------------------------- 
    520558  !- 
     559  !- Some default values so that the operations before the ORCHIDEE initialisation do not fail. 
     560  !- 
     561  z0m(:) = 0.1 
     562  albedo(:,:) = 0.13 
     563  !- 
    521564  itau = 0 
    522   julian = date0 + (itau+0.5)*(dt/one_day) 
    523   CALL ju2ymds (julian, year, month, day, sec) 
    524565  !- 
    525566  !--------------------------------------------------------------------------------------- 
    … …  
    532573  DO itau=1,nbdt 
    533574     ! 
    534      julian = date0 + (itau-0.5)*(dt/one_day) 
    535      ! Needed for STOMATE but should be taken from the arguments of SECHIBA 
    536      in_julian = itau2date(itau, date0, dt) 
    537      ! 
    538      CALL ju2ymds (julian, year, month, day, sec) 
    539      CALL ymds2ju (year,1,1,zero, julian0) 
    540      julian_diff = in_julian-julian0 
     575     CALL time_nextstep( itau ) 
     576     julian = (julian_start + julian_end) /2.0  !julian_end 
    541577     ! 
    542578     ! 
    … …  
    581617     peqAcoef(:) = zero 
    582618     ! 
    583      CALL orchoasis_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, z0, & 
     619     CALL orchoasis_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, z0m, & 
    584620          &                    "Z0 before compute", ktest) 
    585621     ! 
    586      ! Interpolate the wind (which is at hight zlev_uv) to the same height  
    587      ! as the temperature and humidity (at zlev_tq). 
    588      ! 
    589      u_tq(:) = u(:)*LOG(zlev_tq(:)/z0(:))/LOG(zlev_uv(:)/z0(:)) 
    590      v_tq(:) = v(:)*LOG(zlev_tq(:)/z0(:))/LOG(zlev_uv(:)/z0(:)) 
    591622     ! 
    592623     CALL orchoasis_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, u_tq, "USED East-ward wind") 
    … …  
    613644     CALL xios_orchidee_update_calendar(itau_sechiba) 
    614645     ! 
     646     ! Use the obtained albedo to diagnose the Downward Solar and the wind at the right level 
     647     ! 
     648     swdown(:) = swnet(:)/(1.-(albedo(:,1)+albedo(:,2))/2.) 
     649     u_tq(:) = u(:)*LOG(zlev_tq(:)/z0m(:))/LOG(zlev_uv(:)/z0m(:)) 
     650     v_tq(:) = v(:)*LOG(zlev_tq(:)/z0m(:))/LOG(zlev_uv(:)/z0m(:)) 
     651     ! 
    615652     IF ( itau == 1 ) THEN 
    616653        ! 
    … …  
    623660        ! 
    624661        CALL sechiba_initialize( & 
    625              itau_sechiba, iim*jjm,      kjpindex,      kindex,                    & 
     662             itau_sechiba, iim*jjm,      kjpindex,      kindex,                     & 
    626663             lalo_loc,     contfrac,     neighbours,    resolution,  zlev_tq,       & 
    627              u_tq,         v_tq,         qair,          temp_air,                  & 
    628              petAcoef,     peqAcoef,     petBcoef,      peqBcoef,                  & 
    629              precip_rain,  precip_snow,  lwdown,        swnet,       swdown,       & 
     664             u_tq,         v_tq,         qair,          temp_air,                   & 
     665             petAcoef,     peqAcoef,     petBcoef,      peqBcoef,                   & 
     666             precip_rain,  precip_snow,  lwdown,        swnet,       swdown,        & 
    630667             pb,           rest_id,      hist_id,       hist2_id,                   & 
    631668             rest_id_stom, hist_id_stom, hist_id_stom_IPCC,                         & 
    632              coastalflow,  riverflow,    tsol_rad,      vevapp,      qsurf,        & 
    633              z0,           albedo,       fluxsens,      fluxlat,     emis,        & 
     669             coastalflow,  riverflow,    tsol_rad,      vevapp,      qsurf,         & 
     670             z0m,          z0h,          albedo,        fluxsens,    fluxlat, emis, & 
    634671             temp_sol_new, cdrag) 
    635672        CALL orchoasis_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, temp_sol_new, "Init temp_sol_new") 
    636         ! 
    637         ! Use the obtained albedo to diagnose the Downward Solar 
    638         ! 
    639         swdown(:) = swnet(:)/(1.-(albedo(:,1)+albedo(:,2))/2.) 
    640673        ! 
    641674        lrestart_read = .FALSE. 
    … …  
    674707          & netco2, carblu, carbwh, carbha, & 
    675708          ! Surface temperatures and surface properties 
    676           & tsol_rad, temp_sol_new, qsurf, albedo, emis, z0, & 
     709          & tsol_rad, temp_sol_new, qsurf, albedo, emis, z0m, z0h, & 
    677710          ! Vegeation, lai and vegetation height 
    678711          & veget, lai, height, & 
    679712          ! File ids 
    680713          & rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC) 
    681      ! 
    682      ! Use the obtained albedo to diagnose the Downward Solar 
    683      ! 
    684      swdown(:) = swnet(:)/(1.-(albedo(:,1)+albedo(:,2))/2.) 
     714!!$     ! 
     715!!$     ! Use the obtained albedo to diagnose the Downward Solar 
     716!!$     ! 
     717!!$     swdown(:) = swnet(:)/(1.-(albedo(:,1)+albedo(:,2))/2.) 
    685718     ! 
    686719     CALL orchoasis_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, temp_sol_new, "Produced temp_sol_new") 
    … …  
    701734     CALL orchoasis_putvar(itau, dt, kjpindex, kindex - (ii_begin - 1),& 
    702735          &                vevapp, fluxsens, fluxlat, coastalflow, riverflow, & 
    703           &                netco2, carblu, tsol_rad, temp_sol_new, qsurf, albedo, emis, z0) 
     736          &                netco2, carblu, tsol_rad, temp_sol_new, qsurf, albedo, emis, z0m, z0h) 
    704737     ! 
    705738     IF ( timemeasure ) THEN 
    … …  
    723756     ! 
    724757     CALL xios_orchidee_send_field("LandPoints" ,(/ ( REAL(ik), ik=1,kjpindex ) /)) 
    725      CALL xios_orchidee_send_field("Areas", area) 
    726      CALL xios_orchidee_send_field("Contfrac",contfrac) 
    727      CALL xios_orchidee_send_field("evap",vevapp*one_day/dt_sechiba) 
    728      CALL xios_orchidee_send_field("evap_alma",vevapp/dt_sechiba) 
    729      CALL xios_orchidee_send_field("coastalflow",coastalflow/dt_sechiba) 
    730      CALL xios_orchidee_send_field("riverflow",riverflow/dt_sechiba) 
    731      CALL xios_orchidee_send_field("temp_sol_C",temp_sol_new-ZeroCelsius) 
    732      CALL xios_orchidee_send_field("temp_sol_K",temp_sol_new) 
    733      CALL xios_orchidee_send_field("fluxsens",fluxsens) 
    734      CALL xios_orchidee_send_field("fluxlat",fluxlat) 
    735      CALL xios_orchidee_send_field("tair",temp_air) 
     758     CALL xios_orchidee_send_field("areas", area) 
     759     CALL xios_orchidee_send_field("contfrac",contfrac) 
     760     CALL xios_orchidee_send_field("temp_air",temp_air) 
    736761     CALL xios_orchidee_send_field("qair",qair) 
    737762     CALL xios_orchidee_send_field("swnet",swnet) 
    738763     CALL xios_orchidee_send_field("swdown",swdown) 
    739764     ! zpb in hPa, output in Pa 
    740      CALL xios_orchidee_send_field("Psurf",pb*100.) 
     765     CALL xios_orchidee_send_field("pb",pb) 
    741766     ! 
    742767     IF ( .NOT. almaoutput ) THEN 
    … …  
    829854  ! 
    830855  CALL histclo 
     856  ! 
    831857  IF(is_root_prc) THEN 
    832858     CALL restclo 
    … …  
    840866  !- 
    841867  CALL orchideeoasis_clear() 
    842   ! 
    843   WRITE(numout,*) "MPI finalized" 
    844868  !- 
    845869  IF ( timemeasure ) THEN 
    … …  
    855879     CALL stop_timer(timer_mpi) 
    856880  ENDIF 
     881  ! 
     882  ! Finalize MPI and XIOS 
     883!!  CALL Finalize_mpi 
     884  CALL xios_context_finalize() 
     885  CALL xios_finalize() 
    857886  !- 
    858887  CALL oasis_terminate(ierror) 
     888!!  CALL MPI_FINALIZE(ierr) 
    859889  ! 
    860890  WRITE(numout,*) "OASIS terminated" 
    … …  
    875905 
    876906    !- Deallocate all variables existing on all procs (list still incomplete) 
    877  
     907     
    878908    IF ( ALLOCATED(lon_glo) ) DEALLOCATE(lon_glo) 
    879909    IF ( ALLOCATED(lat_glo) ) DEALLOCATE(lat_glo) 
    … …  
    925955    IF ( ALLOCATED(albedo) ) DEALLOCATE(albedo) 
    926956    IF ( ALLOCATED(emis) ) DEALLOCATE(emis) 
    927     IF ( ALLOCATED(z0) ) DEALLOCATE(z0) 
     957 
     958!!$    IF ( ALLOCATED(z0m) ) DEALLOCATE(z0m) 
     959!!$    IF ( ALLOCATED(z0h) ) DEALLOCATE(z0h) 
    928960 
    929961    CALL sechiba_clear() 
Note: See TracChangeset for help on using the changeset viewer.