! ! $Id: interp_sulf_input.f90 5559 2025-02-24 14:39:32Z lebasn $ ! SUBROUTINE interp_sulf_input(debutphy,pdtphys,paprs,tr_seri) USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_dimid, & nf95_inq_varid, nf95_inquire_dimension, nf95_open USE netcdf, ONLY: nf90_get_var, nf90_noerr, nf90_nowrite USE mod_grid_phy_lmdz USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root USE mod_phys_lmdz_omp_data, ONLY : is_omp_root USE phys_local_var_mod, ONLY : budg_3D_backgr_ocs, budg_3D_backgr_so2 USE phys_local_var_mod, ONLY : OCS_lifetime, SO2_lifetime, H2SO4_lifetime, O3_clim USE mod_phys_lmdz_para USE dimphy USE phys_cal_mod USE infotrac_phy USE aerophys USE yomcst_mod_h USE strataer_local_var_mod, ONLY : flag_newclim_file,flag_verbose_strataer IMPLICIT NONE ! Variable input REAL paprs(klon,klev+1) REAL tr_seri(klon,klev,nbtr) REAL, INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde) LOGICAL, INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique ! Variables locales INTEGER n_lat ! number of latitudes in the input data INTEGER n_lon ! number of longitudes in the input data INTEGER, SAVE :: n_lev ! number of levels in the input data !$OMP THREADPRIVATE(n_lev) INTEGER n_mth ! number of months in the input data INTEGER, SAVE :: mth_pre !$OMP THREADPRIVATE(mth_pre) ! Champs reconstitues REAL paprs_glo(klon_glo,klev+1) REAL, ALLOCATABLE :: latitude(:) ! (of input data sorted in strictly ascending order) REAL, ALLOCATABLE :: longitude(:) ! (of input data sorted in strictly ascending order) REAL, ALLOCATABLE :: time(:) ! (of input data sorted in strictly ascending order) REAL, ALLOCATABLE :: lev(:) ! levels of input data REAL, ALLOCATABLE :: OCS_clim_in(:, :, :, :) REAL, ALLOCATABLE :: SO2_clim_in(:, :, :, :) REAL, ALLOCATABLE :: OCS_clim_mth(:, :, :) REAL, ALLOCATABLE :: SO2_clim_mth(:, :, :) REAL, ALLOCATABLE :: OCS_clim_tmp(:, :) REAL, ALLOCATABLE :: SO2_clim_tmp(:, :) REAL OCS_clim_glo(klon_glo,klev) REAL SO2_clim_glo(klon_glo,klev) REAL, ALLOCATABLE :: O3_clim_in(:, :, :, :) REAL, ALLOCATABLE :: O3_clim_mth(:, :, :) REAL, ALLOCATABLE :: O3_clim_tmp(:, :) REAL O3_clim_glo(klon_glo,klev) ! fixed climatos REAL, ALLOCATABLE :: OCS_lifetime_in(:, :, :, :) REAL, ALLOCATABLE :: SO2_lifetime_in(:, :, :, :) REAL, ALLOCATABLE :: OCS_lifetime_mth(:, :, :) REAL, ALLOCATABLE :: SO2_lifetime_mth(:, :, :) REAL, ALLOCATABLE :: OCS_lifetime_tmp(:, :) REAL, ALLOCATABLE :: SO2_lifetime_tmp(:, :) REAL OCS_lifetime_glo(klon_glo,klev) REAL SO2_lifetime_glo(klon_glo,klev) REAL, ALLOCATABLE :: H2SO4_lifetime_in(:, :, :, :) REAL, ALLOCATABLE :: H2SO4_lifetime_mth(:, :, :) REAL, ALLOCATABLE :: H2SO4_lifetime_tmp(:, :) REAL H2SO4_lifetime_glo(klon_glo,klev) ! REAL, ALLOCATABLE, SAVE :: OCS_clim(:,:) REAL, ALLOCATABLE, SAVE :: SO2_clim(:,:) !$OMP THREADPRIVATE(OCS_clim,SO2_clim) ! INTEGER i, k, kk, j REAL p_bound ! For NetCDF: INTEGER :: ncid_in ! IDs for input files INTEGER :: varid, ncerr CHARACTER (LEN=3) :: nc_lat, nc_lon CHARACTER (LEN=256) :: nc_fname INTEGER, PARAMETER :: lev_input=17 !--pressure at interfaces of input data (in Pa) REAL, DIMENSION(lev_input+1), PARAMETER :: & paprs_input=(/ & 1.00000002e+05, 6.06530673e+04, 3.67879449e+04, & 2.23130165e+04, 1.35335286e+04, 8.20850004e+03, & 4.97870695e+03, 3.01973841e+03, 1.83156393e+03, & 1.11089968e+03, 6.73794715e+02, 4.08677153e+02, & 2.47875223e+02, 1.50343923e+02, 9.11881985e+01, & 5.53084382e+01, 3.35462635e+01, 0.0 /) ! IF (.NOT.ALLOCATED(OCS_clim)) ALLOCATE(OCS_clim(klon,klev)) IF (.NOT.ALLOCATED(SO2_clim)) ALLOCATE(SO2_clim(klon,klev)) IF (debutphy.OR.mth_cur.NE.mth_pre) THEN !--preparation of global fields CALL gather(paprs, paprs_glo) IF (is_mpi_root.AND.is_omp_root) THEN OCS_lifetime(:,:)=0.0 SO2_lifetime(:,:)=0.0 H2SO4_lifetime(:,:)=0.0 O3_clim(:,:)=0.0 !--init ncdf variables IF(flag_newclim_file) THEN nc_fname = "ocs_so2_h2so4_annual_lmdz.nc" nc_lat = "LAT" nc_lon = "LON" ELSE ! old file for retro compatibility nc_fname = "ocs_so2_annual_lmdz.nc" nc_lat = "lat" nc_lon = "lon" ENDIF !--reading emission files CALL nf95_open(nc_fname, nf90_nowrite, ncid_in) IF(flag_verbose_strataer) print *, "Open file ", nc_fname CALL nf95_inq_varid(ncid_in, "LEV", varid) CALL nf95_gw_var(ncid_in, varid, lev) n_lev = size(lev) CALL nf95_inq_varid(ncid_in, nc_lat, varid) CALL nf95_gw_var(ncid_in, varid, latitude) n_lat = size(latitude) CALL nf95_inq_varid(ncid_in, nc_lon, varid) CALL nf95_gw_var(ncid_in, varid, longitude) n_lon = size(longitude) CALL nf95_inq_varid(ncid_in, "TIME", varid) CALL nf95_gw_var(ncid_in, varid, time) n_mth = size(time) IF (.NOT.ALLOCATED(OCS_clim_in)) ALLOCATE(OCS_clim_in(n_lon, n_lat, n_lev, n_mth)) IF (.NOT.ALLOCATED(SO2_clim_in)) ALLOCATE(SO2_clim_in(n_lon, n_lat, n_lev, n_mth)) IF (.NOT.ALLOCATED(O3_clim_in)) ALLOCATE(O3_clim_in(n_lon, n_lat, n_lev, n_mth)) IF (.NOT.ALLOCATED(OCS_lifetime_in)) ALLOCATE(OCS_lifetime_in(n_lon, n_lat, n_lev, n_mth)) IF (.NOT.ALLOCATED(SO2_lifetime_in)) ALLOCATE(SO2_lifetime_in(n_lon, n_lat, n_lev, n_mth)) IF (.NOT.ALLOCATED(H2SO4_lifetime_in)) ALLOCATE(H2SO4_lifetime_in(n_lon, n_lat, n_lev, n_mth)) CALL nf95_inq_varid(ncid_in, "OCS", varid) ncerr = nf90_get_var(ncid_in, varid, OCS_clim_in) IF(flag_verbose_strataer) print *,'code erreur OCS=', ncerr, varid CALL nf95_inq_varid(ncid_in, "SO2", varid) ncerr = nf90_get_var(ncid_in, varid, SO2_clim_in) IF(flag_verbose_strataer) print *,'code erreur SO2=', ncerr, varid CALL nf95_inq_varid(ncid_in, "OCS_LIFET", varid) ncerr = nf90_get_var(ncid_in, varid, OCS_lifetime_in) IF(flag_verbose_strataer) print *,'code erreur OCS_lifetime_in=', ncerr, varid CALL nf95_inq_varid(ncid_in, "SO2_LIFET", varid) ncerr = nf90_get_var(ncid_in, varid, SO2_lifetime_in) IF(flag_verbose_strataer) print *,'code erreur SO2_lifetime_in=', ncerr, varid IF(flag_newclim_file) THEN CALL nf95_inq_varid(ncid_in, "O3", varid) ncerr = nf90_get_var(ncid_in, varid, O3_clim_in) IF(flag_verbose_strataer) print *,'code erreur O3=', ncerr, varid CALL nf95_inq_varid(ncid_in, "H2SO4_LIFET", varid) ncerr = nf90_get_var(ncid_in, varid, H2SO4_lifetime_in) IF(flag_verbose_strataer) print *,'code erreur H2SO4_lifetime_in=', ncerr, varid ENDIF CALL nf95_close(ncid_in) IF (.NOT.ALLOCATED(OCS_clim_mth)) ALLOCATE(OCS_clim_mth(n_lon, n_lat, n_lev)) IF (.NOT.ALLOCATED(SO2_clim_mth)) ALLOCATE(SO2_clim_mth(n_lon, n_lat, n_lev)) IF (.NOT.ALLOCATED(OCS_clim_tmp)) ALLOCATE(OCS_clim_tmp(klon_glo, n_lev)) IF (.NOT.ALLOCATED(SO2_clim_tmp)) ALLOCATE(SO2_clim_tmp(klon_glo, n_lev)) IF (.NOT.ALLOCATED(O3_clim_mth)) ALLOCATE(O3_clim_mth(n_lon, n_lat, n_lev)) IF (.NOT.ALLOCATED(O3_clim_tmp)) ALLOCATE(O3_clim_tmp(klon_glo, n_lev)) IF (.NOT.ALLOCATED(OCS_lifetime_mth)) ALLOCATE(OCS_lifetime_mth(n_lon, n_lat, n_lev)) IF (.NOT.ALLOCATED(SO2_lifetime_mth)) ALLOCATE(SO2_lifetime_mth(n_lon, n_lat, n_lev)) IF (.NOT.ALLOCATED(OCS_lifetime_tmp)) ALLOCATE(OCS_lifetime_tmp(klon_glo, n_lev)) IF (.NOT.ALLOCATED(SO2_lifetime_tmp)) ALLOCATE(SO2_lifetime_tmp(klon_glo, n_lev)) IF (.NOT.ALLOCATED(H2SO4_lifetime_mth)) ALLOCATE(H2SO4_lifetime_mth(n_lon, n_lat, n_lev)) IF (.NOT.ALLOCATED(H2SO4_lifetime_tmp)) ALLOCATE(H2SO4_lifetime_tmp(klon_glo, n_lev)) !---select the correct month, undo multiplication with 1.e12 (precision reasons) !---correct latitudinal order and convert input from volume mixing ratio to mass mixing ratio DO j=1,n_lat SO2_clim_mth(:,j,:) = 1.e-12*SO2_clim_in(:,n_lat+1-j,:,mth_cur)*mSO2mol/mAIRmol OCS_clim_mth(:,j,:) = 1.e-12*OCS_clim_in(:,n_lat+1-j,:,mth_cur)*mOCSmol/mAIRmol SO2_lifetime_mth(:,j,:) = SO2_lifetime_in(:,n_lat+1-j,:,mth_cur) OCS_lifetime_mth(:,j,:) = OCS_lifetime_in(:,n_lat+1-j,:,mth_cur) ! O3 from 2d model is not tracer, in VMR IF(flag_newclim_file) THEN H2SO4_lifetime_mth(:,j,:) = H2SO4_lifetime_in(:,n_lat+1-j,:,mth_cur) ! new input files O3_clim_mth(:,j,:) = 1.e-6*O3_clim_in(:,n_lat+1-j,:,mth_cur) ELSE H2SO4_lifetime_mth(:,j,:) = 1.e-6 O3_clim_mth(:,j,:) = 1.e-6 ENDIF ENDDO !---reduce to a klon_glo grid but keep the levels CALL grid2dTo1d_glo(OCS_clim_mth,OCS_clim_tmp) CALL grid2dTo1d_glo(SO2_clim_mth,SO2_clim_tmp) CALL grid2dTo1d_glo(O3_clim_mth,O3_clim_tmp) CALL grid2dTo1d_glo(OCS_lifetime_mth,OCS_lifetime_tmp) CALL grid2dTo1d_glo(SO2_lifetime_mth,SO2_lifetime_tmp) CALL grid2dTo1d_glo(H2SO4_lifetime_mth,H2SO4_lifetime_tmp) !--set lifetime to very high value in uninsolated areas DO i=1, klon_glo DO kk=1, n_lev !--added tests on VMR of species for realism IF (OCS_clim_tmp(i,kk).LT.1.e-20) THEN OCS_clim_tmp(i,kk)=1.0e-20 ENDIF IF (SO2_clim_tmp(i,kk).LT.1.e-20) THEN SO2_clim_tmp(i,kk)=1.0e-20 ENDIF ! 50 ppbv min for O3 IF (O3_clim_tmp(i,kk).LT.50.0e-9) THEN O3_clim_tmp(i,kk)=50.0e-9 ENDIF IF (OCS_lifetime_tmp(i,kk)==0.0) THEN OCS_lifetime_tmp(i,kk)=1.0e12 ENDIF IF (SO2_lifetime_tmp(i,kk)==0.0) THEN SO2_lifetime_tmp(i,kk)=1.0e12 ENDIF IF (H2SO4_lifetime_tmp(i,kk)==0.0) THEN H2SO4_lifetime_tmp(i,kk)=1.0e12 ENDIF ENDDO ENDDO !---regrid weighted lifetime and climatologies DO i=1, klon_glo DO k=1, klev OCS_lifetime_glo(i,k)=0.0 SO2_lifetime_glo(i,k)=0.0 O3_clim_glo(i,k)=0.0 OCS_clim_glo(i,k)=0.0 SO2_clim_glo(i,k)=0.0 H2SO4_lifetime_glo(i,k)=0.0 DO kk=1, n_lev OCS_lifetime_glo(i,k)=OCS_lifetime_glo(i,k)+ & MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) & *OCS_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1)) SO2_lifetime_glo(i,k)=SO2_lifetime_glo(i,k)+ & MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) & *SO2_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1)) IF(flag_newclim_file) THEN H2SO4_lifetime_glo(i,k)=H2SO4_lifetime_glo(i,k)+ & MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk)) & -MAX(paprs_glo(i,k+1),paprs_input(kk+1))) & *H2SO4_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1)) ENDIF OCS_clim_glo(i,k)=OCS_clim_glo(i,k)+ & MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) & *OCS_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1)) SO2_clim_glo(i,k)=SO2_clim_glo(i,k)+ & MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) & *SO2_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1)) IF(flag_newclim_file) THEN O3_clim_glo(i,k)=O3_clim_glo(i,k)+ & MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk)) & -MAX(paprs_glo(i,k+1),paprs_input(kk+1))) & *O3_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1)) ENDIF ENDDO ENDDO ENDDO ENDIF !--is_mpi_root and is_omp_root !--keep memory of previous month mth_pre=mth_cur !--scatter global fields around CALL scatter(OCS_clim_glo, OCS_clim) CALL scatter(SO2_clim_glo, SO2_clim) CALL scatter(O3_clim_glo, O3_clim) CALL scatter(OCS_lifetime_glo, OCS_lifetime) CALL scatter(SO2_lifetime_glo, SO2_lifetime) CALL scatter(H2SO4_lifetime_glo, H2SO4_lifetime) IF (is_mpi_root.AND.is_omp_root) THEN ! DEALLOCATE(OCS_clim_in,SO2_clim_in,O3_clim_in) DEALLOCATE(OCS_clim_mth,SO2_clim_mth,O3_clim_mth) DEALLOCATE(OCS_clim_tmp,SO2_clim_tmp,O3_clim_tmp) DEALLOCATE(OCS_lifetime_in, SO2_lifetime_in, H2SO4_lifetime_in) DEALLOCATE(OCS_lifetime_mth, SO2_lifetime_mth, H2SO4_lifetime_mth) DEALLOCATE(OCS_lifetime_tmp, SO2_lifetime_tmp, H2SO4_lifetime_tmp) ! ENDIF !--is_mpi_root and is_omp_root ENDIF ! debutphy.OR.new month !--set to background value everywhere in the very beginning, later only in the troposphere !--a little dangerous as the MAXVAL is not computed on the global field IF (debutphy.AND.MAXVAL(tr_seri).LT.1.e-30) THEN p_bound=0.0 ELSE p_bound=50000. ENDIF !--regridding tracer concentration on the vertical budg_3D_backgr_ocs(:,:)=0.0 budg_3D_backgr_so2(:,:)=0.0 DO i=1, klon DO k=1, klev ! !--OCS and SO2 prescribed back to their clim values below p_bound IF (paprs(i,k).GT.p_bound) THEN budg_3D_backgr_ocs(i,k)=OCS_clim(i,k)-tr_seri(i,k,id_OCS_strat) budg_3D_backgr_so2(i,k)=SO2_clim(i,k)-tr_seri(i,k,id_SO2_strat) tr_seri(i,k,id_OCS_strat)=OCS_clim(i,k) tr_seri(i,k,id_SO2_strat)=SO2_clim(i,k) ENDIF ENDDO ENDDO !convert SO2_backgr_tend from kg(SO2)/kgA to kg(S)/m2/layer/s for saving as diagnostic DO i=1, klon DO k=1, klev budg_3D_backgr_ocs(i,k)=budg_3D_backgr_ocs(i,k)*mSatom/mOCSmol*(paprs(i,k)-paprs(i,k+1))/RG/pdtphys budg_3D_backgr_so2(i,k)=budg_3D_backgr_so2(i,k)*mSatom/mSO2mol*(paprs(i,k)-paprs(i,k+1))/RG/pdtphys ENDDO ENDDO RETURN END SUBROUTINE interp_sulf_input