PROGRAM main_bpc_hpg !-------------------------------------------------------------------------------------------- ! 1. Declarations USE par_kind USE netcdf USE len_oce USE phycst USE in_out_manager ! I/O manager USE eosinsitu USE zpshde USE dynhpgs USE zdfmxl USE ldfslp USE traldf_iso USE traldf_iso_tile USE tiletype IMPLICIT NONE INTEGER :: ji, jj, jk, jt INTEGER :: jpi_in, jpj_in, jpk_in INTEGER :: istep,maxstep INTEGER :: iostat, ncid, dimid, varid, chunksize INTEGER :: x_dimid, y_dimid, z_dimid, i_dimid, h_dimid INTEGER :: x_varid, y_varid, z_varid, nav_lon_varid, nav_lat_varid, nav_lev_varid INTEGER :: hpgi_varid, hpgj_varid, ta_varid, sa_varid, hbpcgi_varid, hbpcgj_varid INTEGER :: numnam = 12 ! unit number for namelist INTEGER :: kpass = 1 ! for simple harmonic isopycnal diffusion INTEGER :: nlb10 = 7 ! w index just below 10 m depth REAL :: fill = -99999.0 ! note that using REAL(wp) here gave errors in NF90_ENDDEF REAL(wp) :: r1_2 = 0.5_wp REAL(wp) :: r2dt = 7200. ! 1 hour timestep in seconds REAL(wp) :: zcoef, zUfac INTEGER :: jsi, jsj, jsk ! tile loop indices (copied into tile) TYPE(tile_type) :: tile ! Fields read in from files REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: in3d REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: in4d REAL(wp), DIMENSION(:,:), ALLOCATABLE :: nav_lon, nav_lat REAL(wp), DIMENSION(:), ALLOCATABLE :: nav_lev REAL(wp), DIMENSION(:,:), ALLOCATABLE :: e1t, e1u, e1v, e2t, e2u, e2v REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: e3t_0, e3u_0, e3v_0, e3w_0 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmask, umask, vmask, wmask REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ts_n, ts_pcbias REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: e3t_n REAL(wp), DIMENSION(:,:), ALLOCATABLE :: sshn REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: hpgi, hpgj, hbpcgi, hbpcgj ! Fields calculated after input files have been read in REAL(wp), DIMENSION(:,:), ALLOCATABLE :: r1_e1u, r1_e2v, e1e2t, r1_e1e2t, e1e2u, e1e2v, r1_e1e2u, r1_e1e2v, e1_e2v, e2_e1u REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: e3u_n, e3v_n, e3w_n, gdept_n, gdepw_n, gde3w_n, rhd, rhop, rn2b, avt, ahtu, ahtv REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: gtsu, gtsv, gtui, gtvi INTEGER, DIMENSION(:,:), ALLOCATABLE :: mbkt, mbku, mbkv, mikt, miku, mikv, nmln REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: ssmask, gru, grv, risfdep, hmld, hmlp, hmlpt, uslpml, vslpml, wslpiml, wslpjml REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: uslp, vslp, wslpi, wslpj, akz, ah_wslp2, omlmask REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ts_b, ts_a, rab_b REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ua, va REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zftu, zftv ! arrays used for diagnostics REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zri, zrj, zhi, zhj ! NB: 3rd dim=1 to use eos REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zti, ztj ! REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zhpi, zhpj CHARACTER(LEN=256) :: meshfile, modelfile, pcbiasfile, outfile, out_diag LOGICAL :: ln_linssh ! True => linear free surface. e3t does not change with time LOGICAL :: ln_zco, ln_zps, ln_sco ! choice of hpg scheme. Only one logical should be true ! REAL(wp):: rn_avt0, rn_Ud LOGICAL :: ln_isfcav ! set to false LOGICAL :: ln_traldf_msc, ln_traldf_blp, ln_traldf_lap LOGICAL :: l_hst, l_ptr ! logicals for diagnostics INTEGER :: ji_len_3D_tile, jj_len_3D_tile, jk_len_3D_tile LOGICAL :: ln_tiled REAL(wp) :: sum_ua,sum_va,sum_ts_a, step_time REAL(wp) :: et !-------------------------------------------------------------------------------------------- ! 2. User inputs ! They include the namelist nam_pcb_hpg ! the namelist nameos in eos_init ! NAMELIST/nam_bpc_hpg/ meshfile, modelfile, pcbiasfile, outfile, out_diag, ln_linssh, ln_zco, ln_zps, ln_sco, jpi_in, jpj_in, jpk_in, & & rn_avt0, rn_Ud, ji_len_3D_tile, jj_len_3D_tile, jk_len_3D_tile, & & maxstep,ln_tiled write(0,*) 'Starting' meshfile = 'mesh_mask.nc' modelfile = 'restart.nc' pcbiasfile = 'pcbias.nc' outfile = 'hpg_Jan_2019.nc' out_diag = 'out_diag_Jan_2019.txt' ln_linssh = .False. ln_zco = .False. ; ln_zps = .False. ; ln_sco = .False. ji_len_3D_tile = 128 ; jj_len_3D_tile = 64 ; jk_len_3D_tile = 3 ! default values for tile lengths rn_avt0 = 1.2E-5 ! values used in experiment u-bd189 (need to check what this is !!!) rn_ud = 0.011 ! values used in experiment u-bd189 ln_isfcav = .False. ln_traldf_msc = .False.; ln_traldf_blp = .False.; ln_traldf_lap = .True. l_hst = .False. ; l_ptr = .False. nlb10 = 7 ln_tiled=.false. maxstep=10 OPEN( UNIT=numnam, FILE='nam_bpc_hpg', FORM='FORMATTED', STATUS='OLD' ) READ( numnam, nam_bpc_hpg ) CLOSE( numnam ) OPEN( UNIT=6, FILE=out_diag, FORM='FORMATTED', STATUS='UNKNOWN' ) WRITE(6, nam_bpc_hpg ) !-------------------------------------------------------------------------------------------- ! 3. Read in the meshmask data chunksize = 32000000 write(0,*) 'Opening mesh' iostat = NF90_OPEN(TRIM(meshfile), MODE=nf90_nowrite, NCID=ncid, CHUNKSIZE=chunksize) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR opening meshfile: ', TRIM(nf90_strerror(iostat)) iostat = NF90_INQ_DIMID(ncid, 'x', dimid) iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpi) iostat = NF90_INQ_DIMID(ncid, 'y', dimid) iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpj) iostat = NF90_INQ_DIMID(ncid, 'z', dimid) iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpk) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading dimensions in meshfile: ', TRIM(nf90_strerror(iostat)) write(0,*) 'jpi, jpj, jpk = ', jpi, jpj, jpk IF ( jpi /= jpi_in .OR. jpj /= jpj_in .OR. jpk /= jpk_in ) THEN write(0,*) ' FATAL ERROR: the input dimensions in the meshmask file and namelist do not match.' write(0,*) ' Change the nam_pcb_hpg namelist or the input file ' STOP END IF jpim1 = jpi - 1 jpjm1 = jpj - 1 jpkm1 = jpk - 1 ALLOCATE( in3d(jpi,jpj,1), in4d(jpi,jpj,jpk,1) ) ALLOCATE( nav_lon(jpi,jpj), nav_lat(jpi,jpj), nav_lev(jpk), & & e1t(jpi,jpj), e1u(jpi,jpj), e1v(jpi,jpj), e2t(jpi,jpj), e2u(jpi,jpj), e2v(jpi,jpj) ) ALLOCATE( e3t_0(jpi,jpj,jpk), e3u_0(jpi,jpj,jpk), e3v_0(jpi,jpj,jpk), e3w_0(jpi,jpj,jpk), & & tmask(jpi,jpj,jpk), umask(jpi,jpj,jpk), vmask(jpi,jpj,jpk), wmask(jpi,jpj,jpk) ) ALLOCATE( ts_n(jpi,jpj,jpk,jpts), ts_pcbias(jpi,jpj,jpk,jpts) ) ALLOCATE( e3t_n(jpi,jpj,jpk), sshn(jpi,jpj) ) ALLOCATE( hpgi(jpi,jpj,jpk), hpgj(jpi,jpj,jpk), hbpcgi(jpi,jpj,jpk), hbpcgj(jpi,jpj,jpk) ) ALLOCATE( zri(jpi,jpj), zrj(jpi,jpj), zhi(jpi,jpj), zhj(jpi,jpj), zti(jpi,jpj,jpts), ztj(jpi,jpj,jpts), & & zhpi(jpi,jpj,jpk), zhpj(jpi,jpj,jpk)) iostat = NF90_INQ_VARID(ncid, 'nav_lon', varid) iostat = NF90_GET_VAR(ncid, varid, nav_lon(:,:)) iostat = NF90_INQ_VARID(ncid, 'nav_lat', varid) iostat = NF90_GET_VAR(ncid, varid, nav_lat(:,:)) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading nav_lon or nav_lat: ', TRIM(nf90_strerror(iostat)) iostat = NF90_INQ_VARID(ncid, 'nav_lev', varid) iostat = NF90_GET_VAR(ncid, varid, nav_lev(:)) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading nav_lev: ', TRIM(nf90_strerror(iostat)) iostat = NF90_INQ_VARID(ncid, 'e1t', varid) iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) !$OMP PARALLEL WORKSHARE e1t(:,:) = in3d(:,:,1) !$OMP END PARALLEL WORKSHARE iostat = NF90_INQ_VARID(ncid, 'e1u', varid) iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) !$OMP PARALLEL WORKSHARE e1u(:,:) = in3d(:,:,1) !$OMP END PARALLEL WORKSHARE iostat = NF90_INQ_VARID(ncid, 'e1v', varid) iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) !$OMP PARALLEL WORKSHARE e1v(:,:) = in3d(:,:,1) !$OMP END PARALLEL WORKSHARE iostat = NF90_INQ_VARID(ncid, 'e2t', varid) iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) !$OMP PARALLEL WORKSHARE e2t(:,:) = in3d(:,:,1) !$OMP END PARALLEL WORKSHARE iostat = NF90_INQ_VARID(ncid, 'e2u', varid) iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) !$OMP PARALLEL WORKSHARE e2u(:,:) = in3d(:,:,1) !$OMP END PARALLEL WORKSHARE iostat = NF90_INQ_VARID(ncid, 'e2v', varid) iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) !$OMP PARALLEL WORKSHARE e2v(:,:) = in3d(:,:,1) !$OMP END PARALLEL WORKSHARE IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading e1, e2 fields: ', TRIM(nf90_strerror(iostat)) iostat = NF90_INQ_VARID(ncid, 'e3t_0', varid) iostat = NF90_GET_VAR(ncid, varid, e3t_0(:,:,:)) iostat = NF90_INQ_VARID(ncid, 'e3u_0', varid) iostat = NF90_GET_VAR(ncid, varid, e3u_0(:,:,:)) iostat = NF90_INQ_VARID(ncid, 'e3v_0', varid) iostat = NF90_GET_VAR(ncid, varid, e3v_0(:,:,:)) iostat = NF90_INQ_VARID(ncid, 'e3w_0', varid) iostat = NF90_GET_VAR(ncid, varid, e3w_0(:,:,:)) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading e3 fields: ', TRIM(nf90_strerror(iostat)) iostat = NF90_INQ_VARID(ncid, 'tmask', varid) iostat = NF90_GET_VAR(ncid, varid, tmask(:,:,:)) iostat = NF90_INQ_VARID(ncid, 'umask', varid) iostat = NF90_GET_VAR(ncid, varid, umask(:,:,:)) iostat = NF90_INQ_VARID(ncid, 'vmask', varid) iostat = NF90_GET_VAR(ncid, varid, vmask(:,:,:)) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading mask fields: ', TRIM(nf90_strerror(iostat)) iostat = NF90_CLOSE(ncid) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR closing mesh file: ', TRIM(nf90_strerror(iostat)) !-------------------------------------------------------------------------------------------- ! 4. Read in the models fields required write(0,*) 'Opening model fields' iostat = NF90_OPEN(TRIM(modelfile), MODE=nf90_nowrite, NCID=ncid, CHUNKSIZE=chunksize) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR opening main fields file: ', TRIM(nf90_strerror(iostat)) iostat = NF90_INQ_DIMID(ncid, 'x', dimid) iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpi) iostat = NF90_INQ_DIMID(ncid, 'y', dimid) iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpj) iostat = NF90_INQ_DIMID(ncid, 'z', dimid) iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpk) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading main field dimensions: ', TRIM(nf90_strerror(iostat)) write(0,*) 'jpi, jpj, jpk = ', jpi, jpj, jpk IF ( jpi /= jpi_in .OR. jpj /= jpj_in .OR. jpk /= jpk_in ) THEN write(0,*) ' FATAL ERROR: the input dimensions in the model fields file and namelist do not match.' write(0,*) ' Change the nam_pcb_hpg namelist or the input file ' STOP END IF iostat = NF90_INQ_VARID(ncid, 'tn', varid) iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:)) !$OMP PARALLEL WORKSHARE ts_n(:,:,:,1) = in4d(:,:,:,1) !$OMP END PARALLEL WORKSHARE IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading ts_n(:,:,:,1) field: ', TRIM(nf90_strerror(iostat)) iostat = NF90_INQ_VARID(ncid, 'sn', varid) iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:)) !$OMP PARALLEL WORKSHARE ts_n(:,:,:,2) = in4d(:,:,:,1) !$OMP END PARALLEL WORKSHARE IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading ts_n(:,:,:,2) field: ', TRIM(nf90_strerror(iostat)) IF ( .NOT.(ln_linssh) ) THEN iostat = NF90_INQ_VARID(ncid, 'e3t_n', varid) iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:)) !$OMP PARALLEL WORKSHARE e3t_n(:,:,:) = in4d(:,:,:,1) !$OMP END PARALLEL WORKSHARE IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading e3t_n field: ', TRIM(nf90_strerror(iostat)) ELSE !$OMP PARALLEL WORKSHARE e3t_n(:,:,:) = e3t_0(:,:,:) !$OMP END PARALLEL WORKSHARE END IF IF ( ln_sco ) THEN iostat = NF90_INQ_VARID(ncid, 'sshn', varid) iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) !$OMP PARALLEL WORKSHARE sshn(:,:) = in3d(:,:,1) !$OMP END PARALLEL WORKSHARE IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading sshn field: ', TRIM(nf90_strerror(iostat)) ELSE !$OMP PARALLEL WORKSHARE sshn(:,:) = 0._wp !$OMP END PARALLEL WORKSHARE END IF iostat = NF90_CLOSE(ncid) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR closing main fields (file: ', TRIM(nf90_strerror(iostat)) write(0,*) 'Opening pcbias fields' iostat = NF90_OPEN(TRIM(pcbiasfile), MODE=nf90_nowrite, NCID=ncid, CHUNKSIZE=chunksize) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR opening pcbias (file: ', TRIM(nf90_strerror(iostat)) iostat = NF90_INQ_VARID(ncid, 'tbias_p', varid) iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:)) !$OMP PARALLEL WORKSHARE ts_pcbias(:,:,:,1) = in4d(:,:,:,1) !$OMP END PARALLEL WORKSHARE IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading ts_pcbias(:,:,:,1) field: ', TRIM(nf90_strerror(iostat)) iostat = NF90_INQ_VARID(ncid, 'sbias_p', varid) iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:)) !$OMP PARALLEL WORKSHARE ts_pcbias(:,:,:,2) = in4d(:,:,:,1) !$OMP END PARALLEL WORKSHARE IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading ts_pcbias(:,:,:,2) field: ', TRIM(nf90_strerror(iostat)) iostat = NF90_CLOSE(ncid) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR closing pcbias (file: ', TRIM(nf90_strerror(iostat)) !-------------------------------------------------------------------------------------------- ! 5. Open the output file and define the dimensions and output names iostat = NF90_CREATE(TRIM(outfile), NF90_CLOBBER, ncid) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR creating file: ', TRIM(nf90_strerror(iostat)) iostat = NF90_DEF_DIM(ncid, 'x', jpi, x_dimid) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining dimension: ', TRIM(nf90_strerror(iostat)) iostat = NF90_DEF_DIM(ncid, 'y', jpj, y_dimid) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining dimension: ', TRIM(nf90_strerror(iostat)) iostat = NF90_DEF_DIM(ncid, 'z', jpk, z_dimid) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining dimension: ', TRIM(nf90_strerror(iostat)) iostat = NF90_DEF_VAR(ncid, 'nav_lon', NF90_FLOAT, (/ x_dimid,y_dimid /), nav_lon_varid) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining nav_lon: ', TRIM(nf90_strerror(iostat)) iostat = NF90_DEF_VAR(ncid, 'nav_lat', NF90_FLOAT, (/ x_dimid,y_dimid /), nav_lat_varid) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining nav_lon: ', TRIM(nf90_strerror(iostat)) iostat = NF90_DEF_VAR(ncid, 'nav_lev', NF90_FLOAT, (/ z_dimid /), nav_lev_varid) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining nav_lon: ', TRIM(nf90_strerror(iostat)) iostat = NF90_DEF_VAR(ncid, 'hpgi', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), hpgi_varid) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining hpgi: ', TRIM(nf90_strerror(iostat)) iostat = NF90_DEF_VAR(ncid, 'hpgj', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), hpgj_varid) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining hpgj: ', TRIM(nf90_strerror(iostat)) iostat = NF90_DEF_VAR(ncid, 'hbpcgi', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), hbpcgi_varid) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining hbpcgi: ', TRIM(nf90_strerror(iostat)) iostat = NF90_DEF_VAR(ncid, 'hbpcgj', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), hbpcgj_varid) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining hbpcgj: ', TRIM(nf90_strerror(iostat)) iostat = NF90_PUT_ATT(ncid, hpgi_varid, '_FillValue', fill) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining fill hpgi_varid : ', TRIM(nf90_strerror(iostat)) iostat = NF90_PUT_ATT(ncid, hpgj_varid, '_FillValue', fill) iostat = NF90_PUT_ATT(ncid, hbpcgi_varid, '_FillValue', fill) iostat = NF90_PUT_ATT(ncid, hbpcgj_varid, '_FillValue', fill) IF(ln_tiled) THEN !iostat = NF90_DEF_VAR(ncid, 'ta', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), ta_varid) !iostat = NF90_DEF_VAR(ncid, 'sa', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), sa_varid) !iostat = NF90_PUT_ATT(ncid, ta_varid, '_FillValue', fill) !iostat = NF90_PUT_ATT(ncid, sa_varid, '_FillValue', fill) ENDIF IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining fill hbpcgj_varid : ', TRIM(nf90_strerror(iostat)) iostat = NF90_ENDDEF(ncid) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR leaving define mode: ', TRIM(nf90_strerror(iostat)) !set timers hpg_zco_time = 0. hpg_zps_time = 0. hpg_sco_time = 0. zps_hde_time = 0. eos_time = 0. eos2d_time = 0. diffusion_time = 0. step_time = 0. !-------------------------------------------------------------------------------------------- ! 6. Additional allocations for variables calculated in later sections #ifndef MANAGE !$ACC ENTER DATA COPYIN(nav_lev, nav_lon, nav_lat, e1t, e1u, e1v, e2t, e2u, e2v,& !$ACC r1_e1u, r1_e2v, e3t_0, e3w_0, tmask, umask, vmask, wmask, sshn, ts_n, ts_pcbias) !$ACC ENTER DATA CREATE(zti, zhi, zri, ztj, zhj, zrj, zhpi, zhpj, & !$ACC rhd, rhop, gtsu, gtsv, mbku, mbkv, gru, grv, gdepw_n, gde3w_n,& !$ACC ua, va, hpgi, hpgj, hbpcgi, hbpcgj, e3t_n, e3w_n, gdept_n, r1_e1u, r1_e2v) #endif ALLOCATE( r1_e1u(jpi,jpj), r1_e2v(jpi,jpj), e1e2t(jpi,jpj), r1_e1e2t(jpi,jpj), e1e2u(jpi,jpj), & & e1e2v(jpi,jpj), r1_e1e2u(jpi,jpj), r1_e1e2v(jpi,jpj), e1_e2v(jpi,jpj), e2_e1u(jpi,jpj) ) ALLOCATE( e3u_n(jpi,jpj,jpk), e3v_n(jpi,jpj,jpk), e3w_n(jpi,jpj,jpk) ) ALLOCATE( gdept_n(jpi,jpj,jpk), gdepw_n(jpi,jpj,jpk), gde3w_n(jpi,jpj,jpk) ) ALLOCATE( rhd(jpi,jpj,jpk), rhop(jpi,jpj,jpk), rn2b(jpi,jpj,jpk), avt(jpi,jpj,jpk), ahtu(jpi,jpj,jpk), ahtv(jpi,jpj,jpk) ) ALLOCATE( gtsu(jpi,jpj, jpts), gtsv(jpi,jpj, jpts), gtui(jpi,jpj,jpts), gtvi(jpi,jpj,jpts) ) ALLOCATE( mbkt(jpi,jpj), mbku(jpi,jpj), mbkv(jpi,jpj), mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), nmln(jpi,jpj) ) ALLOCATE( ssmask(jpi,jpj), gru(jpi,jpj), grv(jpi,jpj), risfdep(jpi,jpj), & hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), & & uslpml(jpi,jpj), vslpml(jpi,jpj), wslpiml(jpi,jpj), wslpjml(jpi,jpj) ) ALLOCATE( uslp(jpi,jpj,jpk), vslp(jpi,jpj,jpk), wslpi(jpi,jpj,jpk), wslpj(jpi,jpj,jpk) ) ALLOCATE( akz(jpi,jpj,jpk), ah_wslp2(jpi,jpj,jpk), omlmask(jpi,jpj,jpk) ) ALLOCATE( ts_b(jpi,jpj,jpk,jpts), ts_a(jpi,jpj,jpk,jpts), rab_b(jpi,jpj,jpk,jpts) ) ALLOCATE( ua(jpi,jpj,jpk), va(jpi,jpj,jpk) ) IF(ln_tiled) ALLOCATE( zftu(jpi,jpj,jpk,jpts), zftv(jpi,jpj,jpk,jpts) ) ! arrays used for diagnostics ! tiled !-------------------------------------------------------------------------------------------- ! 6. Calculate mbkt, mbku, mbkv, ice-shelf fields, e3t_n, e3w_n, gdept, r1_e1u, r1_e2v !$ACC KERNELS !$OMP PARALLEL !$OMP WORKSHARE mbkt(:,:) = 0 ; mbku(:,:) = 0 ; mbkv(:,:) = 0 !$OMP END WORKSHARE !$OMP DO DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi IF ( tmask(ji,jj,jk) == 1._wp ) mbkt(ji,jj) = jk IF ( umask(ji,jj,jk) == 1._wp ) mbku(ji,jj) = jk IF ( vmask(ji,jj,jk) == 1._wp ) mbkv(ji,jj) = jk END DO END DO END DO !$OMP END DO !$OMP DO DO jj = 1, jpj DO ji = 1, jpi mbkt(ji,jj) = MAX( mbkt(ji,jj) , 1) mbku(ji,jj) = MAX( mbku(ji,jj) , 1) mbkv(ji,jj) = MAX( mbkv(ji,jj) , 1) END DO END DO !$OMP END DO ! set ice-shelf fields for no ice ! mikt(:,:) = 1 ; miku(:,:) = 1 ; mikv(:,:) = 1 ; risfdep(:,:) = 0._wp gtui(:,:,:) = 0._wp ; gtvi(:,:,:) = 0._wp ! The following code is taken from dommsk: dom_msk !$OMP WORKSHARE wmask (:,:,1) = tmask(:,:,1) !$OMP END WORKSHARE !$OMP DO DO jk = 2, jpk ! interior values wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1) END DO !$OMP END DO ! The following code is taken from domvvl: dom_vvl_interpol for case 'W' interpolation IF ( ln_linssh ) THEN !$OMP WORKSHARE e3w_n(:,:,:) = e3w_0(:,:,:) !$OMP END WORKSHARE ELSE !$OMP WORKSHARE e3w_n(:,:,1) = e3w_0(:,:,1) + e3t_n(:,:,1) - e3t_0(:,:,1) !$OMP END WORKSHARE ! - ML - The use of mask in this formula enables the special treatment of the last w-point without indirect adressing !$OMP DO DO jk = 2, jpk e3w_n(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( e3t_n(:,:,jk-1) - e3t_0(:,:,jk-1) ) & & + 0.5_wp * tmask(:,:,jk) * ( e3t_n(:,:,jk ) - e3t_0(:,:,jk ) ) END DO !$OMP END DO END IF ! The following code is taken from domhgr: SUBROUTINE dom_hgr !$OMP WORKSHARE e1e2t (:,:) = e1t(:,:) * e2t(:,:) ; r1_e1e2t(:,:) = 1._wp / e1e2t(:,:) e1e2u (:,:) = e1u(:,:) * e2u(:,:) ; r1_e1e2u(:,:) = 1._wp / e1e2u(:,:) e1e2v (:,:) = e1v(:,:) * e2v(:,:) ; r1_e1e2v(:,:) = 1._wp / e1e2v(:,:) e2_e1u(:,:) = e2u(:,:) / e1u(:,:) ; e1_e2v(:,:) = e1v(:,:) / e2v(:,:) !$OMP END WORKSHARE ! The following code is taken from domvvl !$OMP DO DO jk = 1, jpk DO jj = 1, jpj-1 DO ji = 1, jpi-1 !* from T- to U-point : hor. surface weighted mean e3u_n(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e1e2u(ji,jj) & & * ( e1e2t(ji ,jj) * ( e3t_n(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & & + e1e2t(ji+1,jj) * ( e3t_n(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) e3u_n(ji,jj,jk) = e3u_n(ji,jj,jk) + e3u_0(ji,jj,jk) !* from T- to V-point : hor. surface weighted mean e3v_n(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e1e2v(ji,jj) & & * ( e1e2t(ji,jj ) * ( e3t_n(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & & + e1e2t(ji,jj+1) * ( e3t_n(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) e3v_n(ji,jj,jk) = e3v_n(ji,jj,jk) + e3v_0(ji,jj,jk) END DO END DO END DO !$OMP END DO ! there needs to be an lbc_lnk call here for e3u_n and e3v_n !$OMP WORKSHARE e3u_n(jpi,:,:) = e3u_n(2,:,:) ; e3v_n(jpi,:,:) = e3v_n(2,:,:) e3u_n(:,jpj,:) = e3u_n(:,2,:) ; e3v_n(:,jpj,:) = e3v_n(:,2,:) !$OMP END WORKSHARE ! The following code is taken from domvvl: dom_vvl_init !$OMP WORKSHARE gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) ! reference to the ocean surface (used for MLD and light penetration) gdepw_n(:,:,1) = 0.0_wp gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) !$OMP END WORKSHARE !$OMP DO PRIVATE(zcoef) DO jk = 2, jpk ! vertical sum DO jj = 1,jpj DO ji = 1,jpi ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) ! ! 0.5 where jk = mikt zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) & & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) END DO END DO END DO !$OMP END DO !$OMP DO DO jj = 1, jpj DO ji = 1, jpi r1_e1u(ji,jj) = 1._wp / e1u(ji,jj) r1_e2v(ji,jj) = 1._wp / e2v(ji,jj) END DO END DO !$OMP END DO !-------------------------------------------------------------------------------------------- ! 7. Initialise fields for calculation of isopycnal diffusion !$OMP WORKSHARE uslp (:,:,:) = 0._wp ; uslpml (:,:) = 0._wp vslp (:,:,:) = 0._wp ; vslpml (:,:) = 0._wp wslpi(:,:,:) = 0._wp ; wslpiml(:,:) = 0._wp wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp akz(:,:,:) = 0._wp ; ah_wslp2(:,:,:) = 0._wp ! just to set values in the halo regions ts_b(:,:,:,:) = ts_n(:,:,:,:) ! set set ts on before time-step to the same values as the now timestep ! simplified choice of vertical diffusivity based on zdfphy avt (:,:, 1 ) = 0._wp ; avt (:,:,jpk) = 0._wp avt(:,:,2:jpkm1) = rn_avt0 !$OMP END WORKSHARE !$OMP END PARALLEL ! simplified choice of horizontal diffusivity based on ldftra: ldf_tra_init ! assumes Laplacian rather than bi-Laplacian isopycnal diffusion (inn = 1) l CASE nn_aht_ijk_t = 20 zUfac = r1_2 *rn_Ud ! from ldfc1d_c2d:ldf_c2d case TRA with knn = 1 !$OMP DO DO jj = 1, jpj DO ji = 1, jpi ahtu(ji,jj,:) = zUfac * MAX( e1u(ji,jj), e2u(ji,jj) ) ahtv(ji,jj,:) = zUfac * MAX( e1v(ji,jj), e2v(ji,jj) ) END DO END DO !$OMP END DO !$ACC END KERNELS ! taken from dommsk !$ACC KERNELS !ssmask (:,:) = MAXVAL( tmask(:,:,:), DIM=3 ) ! Causes a run time failure on GPU, work around below !$OMP DO DO jj = 1, jpj DO ji = 1, jpi ssmask(ji,jj)=MAXVAL(tmask(ji,jj,:), DIM=1) END DO END DO !$OMP END DO !$ACC END KERNELS !-------------------------------------------------------------------------------------------- ! 8. Calculate the pressure gradients CALL eos_init ! IF istep == 1 calculate pressure gradients without the bias pressure correction (bpc) ; if istep == 2 calculate them with the bpc WRITE(0,*)"Stepping ",maxstep,"steps" IF ( ln_zco ) THEN write(0,*) ' Calling hpg_zco' ELSE IF ( ln_zps ) THEN write(0,*) ' Calling hpg_zps' ELSE IF ( ln_sco ) THEN write(0,*) ' Calling hpg_sco' END IF IF(ln_tiled) THEN WRITE(0,*)" Tiled isopycnal diffusion",ji_len_3D_tile, jj_len_3D_tile, jk_len_3D_tile ELSE WRITE(0,*)" Non-Tiled isopycnal diffusion" ENDIF WRITE(0,*) step_time = TIMER() DO istep = 1, maxstep !$ACC KERNELS IF ( istep == maxstep ) THEN !$OMP PARALLEL WORKSHARE ts_n(:,:,:,:) = ts_n(:,:,:,:) + ts_pcbias(:,:,:,:) !$OMP END PARALLEL WORKSHARE END IF !$OMP PARALLEL WORKSHARE ua(:,:,:) = 0.0_wp ; va(:,:,:) = 0.0_wp !$OMP END PARALLEL WORKSHARE !$ACC END KERNELS et = TIMER() CALL eos_insitu_pot( ts_n, tmask, rhd, rhop, gdept_n ) ! now in situ density for hpg computation; ts, tmask, gdept are IN; rhd, rhop are OUT eos_time = eos_time + (TIMER() - et) IF( ln_zps ) THEN CALL zps_hde( istep, jpts, ts_n, mbku, mbkv, e3w_n, gdept_n, tmask, umask, vmask, & ! Partial steps: before horizontal gradient & gtsu, gtsv, rhd, gru , grv, & ! of t, s, rd at the last ocean level & zti, zhi, zri, ztj, zhj, zrj ) END IF IF ( ln_zco ) THEN ! z-coordinate et = TIMER() CALL hpg_zco( grav, r1_e1u, r1_e2v, e3w_n, rhd, ua, va, zhpi ,zhpj) hpg_zco_time = hpg_zco_time + (TIMER() - et) ELSE IF ( ln_zps ) THEN ! z-coordinate plus partial steps (interpolation) et = TIMER() CALL hpg_zps( grav, r1_e1u, r1_e2v, mbku, mbkv, gru, grv, e3w_n, rhd, ua, va, zhpi ,zhpj) hpg_zps_time = hpg_zps_time + (TIMER() - et) ELSE IF ( ln_sco ) THEN et = TIMER() CALL hpg_sco( grav, r1_e1u, r1_e2v, e3w_n, rhd, gde3w_n, ln_linssh, ua, va, zhpi ,zhpj ) hpg_sco_time = hpg_sco_time + (TIMER() - et) END IF !$ACC KERNELS IF (istep == 1) THEN !$OMP PARALLEL WORKSHARE hpgi(:,:,:) = ua(:,:,:) hpgj(:,:,:) = va(:,:,:) !$OMP END PARALLEL WORKSHARE ELSE !$OMP PARALLEL WORKSHARE hbpcgi(:,:,:) = hpgi(:,:,:) - ua(:,:,:) hbpcgj(:,:,:) = hpgj(:,:,:) - va(:,:,:) !$OMP END PARALLEL WORKSHARE END IF !$OMP PARALLEL DO DO jk = 2, jpk ! vertical sum DO jj = 1,jpj DO ji = 1,jpi IF ( umask(ji,jj,jk) == 0. ) ua(ji,jj,jk) = fill IF ( vmask(ji,jj,jk) == 0. ) va(ji,jj,jk) = fill END DO END DO END DO !$OMP END PARALLEL DO !$ACC END KERNELS !$ACC KERNELS !$OMP PARALLEL WORKSHARE ! This doesn't seem to scale for OMP sum_ua=SUM(ua) sum_va=SUM(va) !$OMP END PARALLEL WORKSHARE !$ACC END KERNELS !-------------------------------------------------------------------------------------------- ! 9. Calculate isopycnal diffusion et = TIMER() CALL eos_rab_3d( ts_b, gdept_n, tmask, rab_b ) !$ACC KERNELS !$OMP PARALLEL WORKSHARE rn2b(:,:,1) = 0._wp ; rn2b(:,:,jpk) = 0._wp !$OMP END PARALLEL WORKSHARE !$ACC END KERNELS CALL bn2 ( ts_b, rab_b, gdepw_n, gdept_n, e3w_n, wmask, rn2b ) CALL zdf_mxl( rn2b, e3w_n, gdept_n, gdepw_n, avt, wmask, ssmask, mbkt, nlb10, nmln, hmld, hmlp, hmlpt ) CALL ldf_slp( ln_zps, ln_isfcav, rhd, rn2b, tmask, umask, vmask, wmask, gru, grv, & & e1t, e2t, r1_e1u, r1_e2v, e3u_n, e3v_n, e3w_n, gdept_n, gdepw_n, mbku, mbkv, & & mikt, miku, mikv, nmln, hmlp, hmlpt, risfdep, & & omlmask, uslpml, vslpml, wslpiml, wslpjml, uslp, vslp, wslpi, wslpj ) IF(.NOT.ln_tiled) THEN CALL tra_ldf_iso( ln_traldf_msc, ln_traldf_blp, ln_traldf_lap, ln_zps, ln_isfcav, jpts, kpass, r2dt, & & e1t, e1u, e1v, e2t, e2u, e2v, e1e2t, e1_e2v, e2_e1u, r1_e1e2t, mbku, mbkv, miku, mikv, & & umask, vmask, wmask, e3t_n, e3u_n, e3v_n, e3w_n, wslpi, wslpj, uslp, vslp, & & ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts_b , ts_b, & & akz, ah_wslp2, ts_a ) ELSE !$OMP PARALLEL DO DO jsk = 1, jpkm1, jk_len_3D_tile tile % jsk = jsk tile % jsk_2 = MAX(jsk, 2) tile % jskm1 = MAX(jsk-1, 1) tile % jek = MIN(jsk + jk_len_3D_tile - 1, jpkm1) tile % jekp1 = MIN(tile % jek + 1, jpkm1) tile % jskh = tile % jsk - nn_hls tile % jekh = tile % jek + nn_hls DO jsj = nn_hls + 1, jpj - nn_hls, jj_len_3D_tile tile % jsj = jsj tile % jsjm1 = jsj-1 tile % jej = MIN(jsj + jj_len_3D_tile - 1, jpj - nn_hls) tile % jejp1 = tile % jej+1 tile % jsjh = tile % jsj - nn_hls tile % jejh = tile % jej + nn_hls DO jsi = nn_hls + 1, jpi - nn_hls, ji_len_3D_tile tile % jsi = jsi tile % jsim1 = jsi-1 tile % jei = MIN(jsi + ji_len_3D_tile - 1, jpi - nn_hls) tile % jeip1 = tile % jei+1 tile % jsih = tile % jsi - nn_hls tile % jeih = tile % jei + nn_hls CALL tra_ldf_iso_tile( tile, ln_traldf_msc, ln_traldf_blp, ln_traldf_lap, ln_zps, ln_isfcav, jpts, kpass, r2dt, & & e1t, e1u, e1v, e2t, e2u, e2v, e1e2t, e1_e2v, e2_e1u, r1_e1e2t, mbku, mbkv, miku, mikv, & & umask, vmask, wmask, e3t_n, e3u_n, e3v_n, e3w_n, wslpi, wslpj, uslp, vslp, & & ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts_b , ts_b, ts_a, & & l_hst, l_ptr, akz, ah_wslp2, zftu, zftv ) END DO ! jis END DO ! jjs END DO ! jks !$OMP END PARALLEL DO ENDIF diffusion_time = diffusion_time + (TIMER() - et) !$ACC KERNELS !$OMP PARALLEL WORKSHARE sum_ts_a=SUM(ts_a(:,:,:,1)) !$OMP END PARALLEL WORKSHARE !$ACC END KERNELS WRITE(0,*) 'istep = ', istep WRITE(0,*) 'UA SUM ',sum_ua WRITE(0,*) 'VA SUM ',sum_va WRITE(0,*) 'TSA SUM ',sum_ts_a WRITE(0,*) 'Step time ', TIMER() - step_time WRITE(0,*) END DO ! istep step_time = TIMER() - step_time !-------------------------------------------------------------------------------------------- ! 10. Write out the pressure gradient fields write(0,*) 'writing fields: after pseudo time-steps; fields are written at tracer lats and lons for simplicity for now ' !iostat = NF90_PUT_VAR(ncid, nav_lon_varid, nav_lon) !iostat = NF90_PUT_VAR(ncid, nav_lat_varid, nav_lat) !iostat = NF90_PUT_VAR(ncid, nav_lev_varid, nav_lev) !iostat = NF90_PUT_VAR(ncid, hpgi_varid, hpgi) !iostat = NF90_PUT_VAR(ncid, hpgj_varid, hpgj) !iostat = NF90_PUT_VAR(ncid, hbpcgi_varid, hbpcgi) !iostat = NF90_PUT_VAR(ncid, hbpcgj_varid, hbpcgj) IF(ln_tiled) THEN ! !iostat = NF90_PUT_VAR(ncid, ta_varid, ts_a(:,:,:,1) ) ! !iostat = NF90_PUT_VAR(ncid, sa_varid, ts_a(:,:,:,2) ) ENDIF !IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR putting variable: ', TRIM(nf90_strerror(iostat)) ! ! Check Validation based on GPU PGI result IF ( ln_zco ) THEN ! z-coordinate WRITE(0,*)"ZCO: Difference in UA",sum_ua+409653003427.6933_wp WRITE(0,*)"ZCO: Difference in VA",sum_va+407945220514.0470_wp ELSE IF ( ln_zps ) THEN ! z-coordinate plus partial steps (interpolation) WRITE(0,*)"ZPS: Difference in UA",sum_ua+409653003427.7428_wp WRITE(0,*)"ZPS: Difference in VA",sum_va+407945220514.0659_wp ELSE IF ( ln_sco ) THEN WRITE(0,*)"SCO: Difference in UA",sum_ua+409653003427.7501_wp WRITE(0,*)"SCO: Difference in VA",sum_va+407945220514.1797_wp END IF #ifndef MANAGE !$ACC UPDATE HOST (hbpcgi, hbpcgj) #endif !-------------------------------------------------------------------------------------------- ! 9. Close the output file and stop ! iostat = NF90_PUT_VAR(ncid, hbpcgi_varid, hbpcgi) iostat = NF90_PUT_VAR(ncid, hbpcgj_varid, hbpcgj) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR putting variable: ', TRIM(nf90_strerror(iostat)) iostat = NF90_CLOSE(ncid) IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR closing file: ', TRIM(nf90_strerror(iostat)) IF ( ln_zco ) THEN ! z-coordinate WRITE(0,*)'hpg_zco CPU time [s] is: ', hpg_zco_time ELSE IF ( ln_zps ) THEN ! z-coordinate plus partial steps (interpolation) WRITE(0,*)'hpg_zps CPU time [s] is: ', hpg_zps_time ELSE IF ( ln_sco ) THEN WRITE(0,*)'hpg_sco CPU time [s] is: ', hpg_sco_time END IF IF( ln_zps ) WRITE(0,*)'zps_hde CPU time [s] is: ', zps_hde_time WRITE(0,*)'eos_time CPU time [s] is: ', eos_time WRITE(0,*)'isopycnal diffusion CPU time [s] is: ', diffusion_time WRITE(0,*)'Total step time CPU time [s] is: ', step_time write(0,*) 'finishing' END PROGRAM main_bpc_hpg