| 1 | PROGRAM main_bpc_hpg |
|---|
| 2 | |
|---|
| 3 | !-------------------------------------------------------------------------------------------- |
|---|
| 4 | ! 1. Declarations |
|---|
| 5 | |
|---|
| 6 | USE par_kind |
|---|
| 7 | USE netcdf |
|---|
| 8 | USE len_oce |
|---|
| 9 | USE phycst |
|---|
| 10 | USE in_out_manager ! I/O manager |
|---|
| 11 | USE eosinsitu |
|---|
| 12 | USE zpshde |
|---|
| 13 | USE dynhpgs |
|---|
| 14 | USE zdfmxl |
|---|
| 15 | USE ldfslp |
|---|
| 16 | USE traldf_iso |
|---|
| 17 | USE traldf_iso_tile |
|---|
| 18 | USE tiletype |
|---|
| 19 | |
|---|
| 20 | IMPLICIT NONE |
|---|
| 21 | |
|---|
| 22 | INTEGER :: ji, jj, jk, jt |
|---|
| 23 | INTEGER :: jpi_in, jpj_in, jpk_in |
|---|
| 24 | INTEGER :: istep,maxstep |
|---|
| 25 | INTEGER :: iostat, ncid, dimid, varid, chunksize |
|---|
| 26 | INTEGER :: x_dimid, y_dimid, z_dimid, i_dimid, h_dimid |
|---|
| 27 | INTEGER :: x_varid, y_varid, z_varid, nav_lon_varid, nav_lat_varid, nav_lev_varid |
|---|
| 28 | INTEGER :: hpgi_varid, hpgj_varid, ta_varid, sa_varid, hbpcgi_varid, hbpcgj_varid |
|---|
| 29 | INTEGER :: numnam = 12 ! unit number for namelist |
|---|
| 30 | INTEGER :: kpass = 1 ! for simple harmonic isopycnal diffusion |
|---|
| 31 | INTEGER :: nlb10 = 7 ! w index just below 10 m depth |
|---|
| 32 | |
|---|
| 33 | REAL :: fill = -99999.0 ! note that using REAL(wp) here gave errors in NF90_ENDDEF |
|---|
| 34 | REAL(wp) :: r1_2 = 0.5_wp |
|---|
| 35 | REAL(wp) :: r2dt = 7200. ! 1 hour timestep in seconds |
|---|
| 36 | REAL(wp) :: zcoef, zUfac |
|---|
| 37 | |
|---|
| 38 | INTEGER :: jsi, jsj, jsk ! tile loop indices (copied into tile) |
|---|
| 39 | TYPE(tile_type) :: tile |
|---|
| 40 | |
|---|
| 41 | ! Fields read in from files |
|---|
| 42 | |
|---|
| 43 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: in3d |
|---|
| 44 | REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: in4d |
|---|
| 45 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: nav_lon, nav_lat |
|---|
| 46 | REAL(wp), DIMENSION(:), ALLOCATABLE :: nav_lev |
|---|
| 47 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: e1t, e1u, e1v, e2t, e2u, e2v |
|---|
| 48 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: e3t_0, e3u_0, e3v_0, e3w_0 |
|---|
| 49 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmask, umask, vmask, wmask |
|---|
| 50 | REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ts_n, ts_pcbias |
|---|
| 51 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: e3t_n |
|---|
| 52 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: sshn |
|---|
| 53 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: hpgi, hpgj, hbpcgi, hbpcgj |
|---|
| 54 | |
|---|
| 55 | ! Fields calculated after input files have been read in |
|---|
| 56 | |
|---|
| 57 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: r1_e1u, r1_e2v, e1e2t, r1_e1e2t, e1e2u, e1e2v, r1_e1e2u, r1_e1e2v, e1_e2v, e2_e1u |
|---|
| 58 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: e3u_n, e3v_n, e3w_n, gdept_n, gdepw_n, gde3w_n, rhd, rhop, rn2b, avt, ahtu, ahtv |
|---|
| 59 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: gtsu, gtsv, gtui, gtvi |
|---|
| 60 | INTEGER, DIMENSION(:,:), ALLOCATABLE :: mbkt, mbku, mbkv, mikt, miku, mikv, nmln |
|---|
| 61 | REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: ssmask, gru, grv, risfdep, hmld, hmlp, hmlpt, uslpml, vslpml, wslpiml, wslpjml |
|---|
| 62 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: uslp, vslp, wslpi, wslpj, akz, ah_wslp2, omlmask |
|---|
| 63 | REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ts_b, ts_a, rab_b |
|---|
| 64 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ua, va |
|---|
| 65 | |
|---|
| 66 | REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zftu, zftv ! arrays used for diagnostics |
|---|
| 67 | REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zri, zrj, zhi, zhj ! NB: 3rd dim=1 to use eos |
|---|
| 68 | REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zti, ztj ! |
|---|
| 69 | REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zhpi, zhpj |
|---|
| 70 | |
|---|
| 71 | CHARACTER(LEN=256) :: meshfile, modelfile, pcbiasfile, outfile, out_diag |
|---|
| 72 | |
|---|
| 73 | LOGICAL :: ln_linssh ! True => linear free surface. e3t does not change with time |
|---|
| 74 | LOGICAL :: ln_zco, ln_zps, ln_sco ! choice of hpg scheme. Only one logical should be true ! |
|---|
| 75 | REAL(wp):: rn_avt0, rn_Ud |
|---|
| 76 | |
|---|
| 77 | LOGICAL :: ln_isfcav ! set to false |
|---|
| 78 | LOGICAL :: ln_traldf_msc, ln_traldf_blp, ln_traldf_lap |
|---|
| 79 | |
|---|
| 80 | LOGICAL :: l_hst, l_ptr ! logicals for diagnostics |
|---|
| 81 | |
|---|
| 82 | INTEGER :: ji_len_3D_tile, jj_len_3D_tile, jk_len_3D_tile |
|---|
| 83 | |
|---|
| 84 | LOGICAL :: ln_tiled |
|---|
| 85 | |
|---|
| 86 | REAL(wp) :: sum_ua,sum_va,sum_ts_a, step_time |
|---|
| 87 | REAL(wp) :: et |
|---|
| 88 | |
|---|
| 89 | !-------------------------------------------------------------------------------------------- |
|---|
| 90 | ! 2. User inputs |
|---|
| 91 | ! They include the namelist nam_pcb_hpg |
|---|
| 92 | ! the namelist nameos in eos_init |
|---|
| 93 | ! |
|---|
| 94 | |
|---|
| 95 | NAMELIST/nam_bpc_hpg/ meshfile, modelfile, pcbiasfile, outfile, out_diag, ln_linssh, ln_zco, ln_zps, ln_sco, jpi_in, jpj_in, jpk_in, & |
|---|
| 96 | & rn_avt0, rn_Ud, ji_len_3D_tile, jj_len_3D_tile, jk_len_3D_tile, & |
|---|
| 97 | & maxstep,ln_tiled |
|---|
| 98 | |
|---|
| 99 | write(0,*) 'Starting' |
|---|
| 100 | |
|---|
| 101 | meshfile = 'mesh_mask.nc' |
|---|
| 102 | modelfile = 'restart.nc' |
|---|
| 103 | pcbiasfile = 'pcbias.nc' |
|---|
| 104 | outfile = 'hpg_Jan_2019.nc' |
|---|
| 105 | out_diag = 'out_diag_Jan_2019.txt' |
|---|
| 106 | ln_linssh = .False. |
|---|
| 107 | ln_zco = .False. ; ln_zps = .False. ; ln_sco = .False. |
|---|
| 108 | |
|---|
| 109 | ji_len_3D_tile = 128 ; jj_len_3D_tile = 64 ; jk_len_3D_tile = 3 ! default values for tile lengths |
|---|
| 110 | |
|---|
| 111 | rn_avt0 = 1.2E-5 ! values used in experiment u-bd189 (need to check what this is !!!) |
|---|
| 112 | rn_ud = 0.011 ! values used in experiment u-bd189 |
|---|
| 113 | |
|---|
| 114 | ln_isfcav = .False. |
|---|
| 115 | |
|---|
| 116 | ln_traldf_msc = .False.; ln_traldf_blp = .False.; ln_traldf_lap = .True. |
|---|
| 117 | |
|---|
| 118 | l_hst = .False. ; l_ptr = .False. |
|---|
| 119 | |
|---|
| 120 | nlb10 = 7 |
|---|
| 121 | |
|---|
| 122 | ln_tiled=.false. |
|---|
| 123 | maxstep=10 |
|---|
| 124 | |
|---|
| 125 | OPEN( UNIT=numnam, FILE='nam_bpc_hpg', FORM='FORMATTED', STATUS='OLD' ) |
|---|
| 126 | READ( numnam, nam_bpc_hpg ) |
|---|
| 127 | CLOSE( numnam ) |
|---|
| 128 | |
|---|
| 129 | OPEN( UNIT=6, FILE=out_diag, FORM='FORMATTED', STATUS='UNKNOWN' ) |
|---|
| 130 | WRITE(6, nam_bpc_hpg ) |
|---|
| 131 | |
|---|
| 132 | !-------------------------------------------------------------------------------------------- |
|---|
| 133 | ! 3. Read in the meshmask data |
|---|
| 134 | |
|---|
| 135 | chunksize = 32000000 |
|---|
| 136 | |
|---|
| 137 | write(0,*) 'Opening mesh' |
|---|
| 138 | |
|---|
| 139 | iostat = NF90_OPEN(TRIM(meshfile), MODE=nf90_nowrite, NCID=ncid, CHUNKSIZE=chunksize) |
|---|
| 140 | |
|---|
| 141 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR opening meshfile: ', TRIM(nf90_strerror(iostat)) |
|---|
| 142 | |
|---|
| 143 | iostat = NF90_INQ_DIMID(ncid, 'x', dimid) |
|---|
| 144 | iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpi) |
|---|
| 145 | iostat = NF90_INQ_DIMID(ncid, 'y', dimid) |
|---|
| 146 | iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpj) |
|---|
| 147 | iostat = NF90_INQ_DIMID(ncid, 'z', dimid) |
|---|
| 148 | iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpk) |
|---|
| 149 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading dimensions in meshfile: ', TRIM(nf90_strerror(iostat)) |
|---|
| 150 | |
|---|
| 151 | write(0,*) 'jpi, jpj, jpk = ', jpi, jpj, jpk |
|---|
| 152 | |
|---|
| 153 | IF ( jpi /= jpi_in .OR. jpj /= jpj_in .OR. jpk /= jpk_in ) THEN |
|---|
| 154 | write(0,*) ' FATAL ERROR: the input dimensions in the meshmask file and namelist do not match.' |
|---|
| 155 | write(0,*) ' Change the nam_pcb_hpg namelist or the input file ' |
|---|
| 156 | STOP |
|---|
| 157 | END IF |
|---|
| 158 | |
|---|
| 159 | jpim1 = jpi - 1 |
|---|
| 160 | jpjm1 = jpj - 1 |
|---|
| 161 | jpkm1 = jpk - 1 |
|---|
| 162 | |
|---|
| 163 | ALLOCATE( in3d(jpi,jpj,1), in4d(jpi,jpj,jpk,1) ) |
|---|
| 164 | |
|---|
| 165 | ALLOCATE( nav_lon(jpi,jpj), nav_lat(jpi,jpj), nav_lev(jpk), & |
|---|
| 166 | & e1t(jpi,jpj), e1u(jpi,jpj), e1v(jpi,jpj), e2t(jpi,jpj), e2u(jpi,jpj), e2v(jpi,jpj) ) |
|---|
| 167 | |
|---|
| 168 | ALLOCATE( e3t_0(jpi,jpj,jpk), e3u_0(jpi,jpj,jpk), e3v_0(jpi,jpj,jpk), e3w_0(jpi,jpj,jpk), & |
|---|
| 169 | & tmask(jpi,jpj,jpk), umask(jpi,jpj,jpk), vmask(jpi,jpj,jpk), wmask(jpi,jpj,jpk) ) |
|---|
| 170 | |
|---|
| 171 | ALLOCATE( ts_n(jpi,jpj,jpk,jpts), ts_pcbias(jpi,jpj,jpk,jpts) ) |
|---|
| 172 | ALLOCATE( e3t_n(jpi,jpj,jpk), sshn(jpi,jpj) ) |
|---|
| 173 | |
|---|
| 174 | ALLOCATE( hpgi(jpi,jpj,jpk), hpgj(jpi,jpj,jpk), hbpcgi(jpi,jpj,jpk), hbpcgj(jpi,jpj,jpk) ) |
|---|
| 175 | |
|---|
| 176 | ALLOCATE( zri(jpi,jpj), zrj(jpi,jpj), zhi(jpi,jpj), zhj(jpi,jpj), zti(jpi,jpj,jpts), ztj(jpi,jpj,jpts), & |
|---|
| 177 | & zhpi(jpi,jpj,jpk), zhpj(jpi,jpj,jpk)) |
|---|
| 178 | |
|---|
| 179 | iostat = NF90_INQ_VARID(ncid, 'nav_lon', varid) |
|---|
| 180 | iostat = NF90_GET_VAR(ncid, varid, nav_lon(:,:)) |
|---|
| 181 | |
|---|
| 182 | iostat = NF90_INQ_VARID(ncid, 'nav_lat', varid) |
|---|
| 183 | iostat = NF90_GET_VAR(ncid, varid, nav_lat(:,:)) |
|---|
| 184 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading nav_lon or nav_lat: ', TRIM(nf90_strerror(iostat)) |
|---|
| 185 | |
|---|
| 186 | iostat = NF90_INQ_VARID(ncid, 'nav_lev', varid) |
|---|
| 187 | iostat = NF90_GET_VAR(ncid, varid, nav_lev(:)) |
|---|
| 188 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading nav_lev: ', TRIM(nf90_strerror(iostat)) |
|---|
| 189 | |
|---|
| 190 | iostat = NF90_INQ_VARID(ncid, 'e1t', varid) |
|---|
| 191 | iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) |
|---|
| 192 | e1t(:,:) = in3d(:,:,1) |
|---|
| 193 | |
|---|
| 194 | iostat = NF90_INQ_VARID(ncid, 'e1u', varid) |
|---|
| 195 | iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) |
|---|
| 196 | e1u(:,:) = in3d(:,:,1) |
|---|
| 197 | |
|---|
| 198 | iostat = NF90_INQ_VARID(ncid, 'e1v', varid) |
|---|
| 199 | iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) |
|---|
| 200 | e1v(:,:) = in3d(:,:,1) |
|---|
| 201 | |
|---|
| 202 | iostat = NF90_INQ_VARID(ncid, 'e2t', varid) |
|---|
| 203 | iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) |
|---|
| 204 | e2t(:,:) = in3d(:,:,1) |
|---|
| 205 | |
|---|
| 206 | iostat = NF90_INQ_VARID(ncid, 'e2u', varid) |
|---|
| 207 | iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) |
|---|
| 208 | e2u(:,:) = in3d(:,:,1) |
|---|
| 209 | |
|---|
| 210 | iostat = NF90_INQ_VARID(ncid, 'e2v', varid) |
|---|
| 211 | iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) |
|---|
| 212 | e2v(:,:) = in3d(:,:,1) |
|---|
| 213 | |
|---|
| 214 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading e1, e2 fields: ', TRIM(nf90_strerror(iostat)) |
|---|
| 215 | |
|---|
| 216 | iostat = NF90_INQ_VARID(ncid, 'e3t_0', varid) |
|---|
| 217 | iostat = NF90_GET_VAR(ncid, varid, e3t_0(:,:,:)) |
|---|
| 218 | |
|---|
| 219 | iostat = NF90_INQ_VARID(ncid, 'e3u_0', varid) |
|---|
| 220 | iostat = NF90_GET_VAR(ncid, varid, e3u_0(:,:,:)) |
|---|
| 221 | |
|---|
| 222 | iostat = NF90_INQ_VARID(ncid, 'e3v_0', varid) |
|---|
| 223 | iostat = NF90_GET_VAR(ncid, varid, e3v_0(:,:,:)) |
|---|
| 224 | |
|---|
| 225 | iostat = NF90_INQ_VARID(ncid, 'e3w_0', varid) |
|---|
| 226 | iostat = NF90_GET_VAR(ncid, varid, e3w_0(:,:,:)) |
|---|
| 227 | |
|---|
| 228 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading e3 fields: ', TRIM(nf90_strerror(iostat)) |
|---|
| 229 | |
|---|
| 230 | iostat = NF90_INQ_VARID(ncid, 'tmask', varid) |
|---|
| 231 | iostat = NF90_GET_VAR(ncid, varid, tmask(:,:,:)) |
|---|
| 232 | |
|---|
| 233 | iostat = NF90_INQ_VARID(ncid, 'umask', varid) |
|---|
| 234 | iostat = NF90_GET_VAR(ncid, varid, umask(:,:,:)) |
|---|
| 235 | |
|---|
| 236 | iostat = NF90_INQ_VARID(ncid, 'vmask', varid) |
|---|
| 237 | iostat = NF90_GET_VAR(ncid, varid, vmask(:,:,:)) |
|---|
| 238 | |
|---|
| 239 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading mask fields: ', TRIM(nf90_strerror(iostat)) |
|---|
| 240 | |
|---|
| 241 | iostat = NF90_CLOSE(ncid) |
|---|
| 242 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR closing mesh file: ', TRIM(nf90_strerror(iostat)) |
|---|
| 243 | |
|---|
| 244 | !-------------------------------------------------------------------------------------------- |
|---|
| 245 | ! 4. Read in the models fields required |
|---|
| 246 | |
|---|
| 247 | write(0,*) 'Opening model fields' |
|---|
| 248 | |
|---|
| 249 | iostat = NF90_OPEN(TRIM(modelfile), MODE=nf90_nowrite, NCID=ncid, CHUNKSIZE=chunksize) |
|---|
| 250 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR opening main fields file: ', TRIM(nf90_strerror(iostat)) |
|---|
| 251 | |
|---|
| 252 | iostat = NF90_INQ_DIMID(ncid, 'x', dimid) |
|---|
| 253 | iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpi) |
|---|
| 254 | iostat = NF90_INQ_DIMID(ncid, 'y', dimid) |
|---|
| 255 | iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpj) |
|---|
| 256 | iostat = NF90_INQ_DIMID(ncid, 'z', dimid) |
|---|
| 257 | iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpk) |
|---|
| 258 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading main field dimensions: ', TRIM(nf90_strerror(iostat)) |
|---|
| 259 | |
|---|
| 260 | write(0,*) 'jpi, jpj, jpk = ', jpi, jpj, jpk |
|---|
| 261 | |
|---|
| 262 | IF ( jpi /= jpi_in .OR. jpj /= jpj_in .OR. jpk /= jpk_in ) THEN |
|---|
| 263 | write(0,*) ' FATAL ERROR: the input dimensions in the model fields file and namelist do not match.' |
|---|
| 264 | write(0,*) ' Change the nam_pcb_hpg namelist or the input file ' |
|---|
| 265 | STOP |
|---|
| 266 | END IF |
|---|
| 267 | |
|---|
| 268 | iostat = NF90_INQ_VARID(ncid, 'tn', varid) |
|---|
| 269 | iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:)) |
|---|
| 270 | ts_n(:,:,:,1) = in4d(:,:,:,1) |
|---|
| 271 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading ts_n(:,:,:,1) field: ', TRIM(nf90_strerror(iostat)) |
|---|
| 272 | |
|---|
| 273 | iostat = NF90_INQ_VARID(ncid, 'sn', varid) |
|---|
| 274 | iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:)) |
|---|
| 275 | ts_n(:,:,:,2) = in4d(:,:,:,1) |
|---|
| 276 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading ts_n(:,:,:,2) field: ', TRIM(nf90_strerror(iostat)) |
|---|
| 277 | |
|---|
| 278 | IF ( .NOT.(ln_linssh) ) THEN |
|---|
| 279 | iostat = NF90_INQ_VARID(ncid, 'e3t_n', varid) |
|---|
| 280 | iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:)) |
|---|
| 281 | e3t_n(:,:,:) = in4d(:,:,:,1) |
|---|
| 282 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading e3t_n field: ', TRIM(nf90_strerror(iostat)) |
|---|
| 283 | ELSE |
|---|
| 284 | e3t_n(:,:,:) = e3t_0(:,:,:) |
|---|
| 285 | END IF |
|---|
| 286 | |
|---|
| 287 | IF ( ln_sco ) THEN |
|---|
| 288 | iostat = NF90_INQ_VARID(ncid, 'sshn', varid) |
|---|
| 289 | iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) |
|---|
| 290 | sshn(:,:) = in3d(:,:,1) |
|---|
| 291 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading sshn field: ', TRIM(nf90_strerror(iostat)) |
|---|
| 292 | ELSE |
|---|
| 293 | sshn(:,:) = 0._wp |
|---|
| 294 | END IF |
|---|
| 295 | |
|---|
| 296 | iostat = NF90_CLOSE(ncid) |
|---|
| 297 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR closing main fields (file: ', TRIM(nf90_strerror(iostat)) |
|---|
| 298 | |
|---|
| 299 | write(0,*) 'Opening pcbias fields' |
|---|
| 300 | |
|---|
| 301 | iostat = NF90_OPEN(TRIM(pcbiasfile), MODE=nf90_nowrite, NCID=ncid, CHUNKSIZE=chunksize) |
|---|
| 302 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR opening pcbias (file: ', TRIM(nf90_strerror(iostat)) |
|---|
| 303 | |
|---|
| 304 | iostat = NF90_INQ_VARID(ncid, 'tbias_p', varid) |
|---|
| 305 | iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:)) |
|---|
| 306 | ts_pcbias(:,:,:,1) = in4d(:,:,:,1) |
|---|
| 307 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading ts_pcbias(:,:,:,1) field: ', TRIM(nf90_strerror(iostat)) |
|---|
| 308 | |
|---|
| 309 | iostat = NF90_INQ_VARID(ncid, 'sbias_p', varid) |
|---|
| 310 | iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:)) |
|---|
| 311 | ts_pcbias(:,:,:,2) = in4d(:,:,:,1) |
|---|
| 312 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading ts_pcbias(:,:,:,2) field: ', TRIM(nf90_strerror(iostat)) |
|---|
| 313 | |
|---|
| 314 | iostat = NF90_CLOSE(ncid) |
|---|
| 315 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR closing pcbias (file: ', TRIM(nf90_strerror(iostat)) |
|---|
| 316 | |
|---|
| 317 | !-------------------------------------------------------------------------------------------- |
|---|
| 318 | ! 5. Open the output file and define the dimensions and output names |
|---|
| 319 | iostat = NF90_CREATE(TRIM(outfile), NF90_CLOBBER, ncid) |
|---|
| 320 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR creating file: ', TRIM(nf90_strerror(iostat)) |
|---|
| 321 | |
|---|
| 322 | iostat = NF90_DEF_DIM(ncid, 'x', jpi, x_dimid) |
|---|
| 323 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining dimension: ', TRIM(nf90_strerror(iostat)) |
|---|
| 324 | iostat = NF90_DEF_DIM(ncid, 'y', jpj, y_dimid) |
|---|
| 325 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining dimension: ', TRIM(nf90_strerror(iostat)) |
|---|
| 326 | iostat = NF90_DEF_DIM(ncid, 'z', jpk, z_dimid) |
|---|
| 327 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining dimension: ', TRIM(nf90_strerror(iostat)) |
|---|
| 328 | |
|---|
| 329 | iostat = NF90_DEF_VAR(ncid, 'nav_lon', NF90_FLOAT, (/ x_dimid,y_dimid /), nav_lon_varid) |
|---|
| 330 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining nav_lon: ', TRIM(nf90_strerror(iostat)) |
|---|
| 331 | iostat = NF90_DEF_VAR(ncid, 'nav_lat', NF90_FLOAT, (/ x_dimid,y_dimid /), nav_lat_varid) |
|---|
| 332 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining nav_lon: ', TRIM(nf90_strerror(iostat)) |
|---|
| 333 | iostat = NF90_DEF_VAR(ncid, 'nav_lev', NF90_FLOAT, (/ z_dimid /), nav_lev_varid) |
|---|
| 334 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining nav_lon: ', TRIM(nf90_strerror(iostat)) |
|---|
| 335 | |
|---|
| 336 | iostat = NF90_DEF_VAR(ncid, 'hpgi', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), hpgi_varid) |
|---|
| 337 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining hpgi: ', TRIM(nf90_strerror(iostat)) |
|---|
| 338 | iostat = NF90_DEF_VAR(ncid, 'hpgj', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), hpgj_varid) |
|---|
| 339 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining hpgj: ', TRIM(nf90_strerror(iostat)) |
|---|
| 340 | iostat = NF90_DEF_VAR(ncid, 'hbpcgi', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), hbpcgi_varid) |
|---|
| 341 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining hbpcgi: ', TRIM(nf90_strerror(iostat)) |
|---|
| 342 | iostat = NF90_DEF_VAR(ncid, 'hbpcgj', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), hbpcgj_varid) |
|---|
| 343 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining hbpcgj: ', TRIM(nf90_strerror(iostat)) |
|---|
| 344 | |
|---|
| 345 | iostat = NF90_PUT_ATT(ncid, hpgi_varid, '_FillValue', fill) |
|---|
| 346 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining fill hpgi_varid : ', TRIM(nf90_strerror(iostat)) |
|---|
| 347 | iostat = NF90_PUT_ATT(ncid, hpgj_varid, '_FillValue', fill) |
|---|
| 348 | iostat = NF90_PUT_ATT(ncid, hbpcgi_varid, '_FillValue', fill) |
|---|
| 349 | iostat = NF90_PUT_ATT(ncid, hbpcgj_varid, '_FillValue', fill) |
|---|
| 350 | IF(ln_tiled) THEN |
|---|
| 351 | !iostat = NF90_DEF_VAR(ncid, 'ta', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), ta_varid) |
|---|
| 352 | !iostat = NF90_DEF_VAR(ncid, 'sa', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), sa_varid) |
|---|
| 353 | !iostat = NF90_PUT_ATT(ncid, ta_varid, '_FillValue', fill) |
|---|
| 354 | !iostat = NF90_PUT_ATT(ncid, sa_varid, '_FillValue', fill) |
|---|
| 355 | ENDIF |
|---|
| 356 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining fill hbpcgj_varid : ', TRIM(nf90_strerror(iostat)) |
|---|
| 357 | |
|---|
| 358 | iostat = NF90_ENDDEF(ncid) |
|---|
| 359 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR leaving define mode: ', TRIM(nf90_strerror(iostat)) |
|---|
| 360 | !set timers |
|---|
| 361 | hpg_zco_time = 0. |
|---|
| 362 | hpg_zps_time = 0. |
|---|
| 363 | hpg_sco_time = 0. |
|---|
| 364 | zps_hde_time = 0. |
|---|
| 365 | eos_time = 0. |
|---|
| 366 | eos2d_time = 0. |
|---|
| 367 | diffusion_time = 0. |
|---|
| 368 | step_time = 0. |
|---|
| 369 | !-------------------------------------------------------------------------------------------- |
|---|
| 370 | ! 6. Additional allocations for variables calculated in later sections |
|---|
| 371 | |
|---|
| 372 | #ifndef MANAGE |
|---|
| 373 | !$ACC ENTER DATA COPYIN(nav_lev, nav_lon, nav_lat, e1t, e1u, e1v, e2t, e2u, e2v,& |
|---|
| 374 | !$ACC r1_e1u, r1_e2v, e3t_0, e3w_0, tmask, umask, vmask, wmask, sshn, ts_n, ts_pcbias) |
|---|
| 375 | |
|---|
| 376 | !$ACC ENTER DATA CREATE(zti, zhi, zri, ztj, zhj, zrj, zhpi, zhpj, & |
|---|
| 377 | !$ACC rhd, rhop, gtsu, gtsv, mbku, mbkv, gru, grv, gdepw_n, gde3w_n,& |
|---|
| 378 | !$ACC ua, va, hpgi, hpgj, hbpcgi, hbpcgj, e3t_n, e3w_n, gdept_n, r1_e1u, r1_e2v) |
|---|
| 379 | #endif |
|---|
| 380 | |
|---|
| 381 | ALLOCATE( r1_e1u(jpi,jpj), r1_e2v(jpi,jpj), e1e2t(jpi,jpj), r1_e1e2t(jpi,jpj), e1e2u(jpi,jpj), & |
|---|
| 382 | & e1e2v(jpi,jpj), r1_e1e2u(jpi,jpj), r1_e1e2v(jpi,jpj), e1_e2v(jpi,jpj), e2_e1u(jpi,jpj) ) |
|---|
| 383 | |
|---|
| 384 | ALLOCATE( e3u_n(jpi,jpj,jpk), e3v_n(jpi,jpj,jpk), e3w_n(jpi,jpj,jpk) ) |
|---|
| 385 | |
|---|
| 386 | ALLOCATE( gdept_n(jpi,jpj,jpk), gdepw_n(jpi,jpj,jpk), gde3w_n(jpi,jpj,jpk) ) |
|---|
| 387 | |
|---|
| 388 | 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) ) |
|---|
| 389 | |
|---|
| 390 | ALLOCATE( gtsu(jpi,jpj, jpts), gtsv(jpi,jpj, jpts), gtui(jpi,jpj,jpts), gtvi(jpi,jpj,jpts) ) |
|---|
| 391 | |
|---|
| 392 | ALLOCATE( mbkt(jpi,jpj), mbku(jpi,jpj), mbkv(jpi,jpj), mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), nmln(jpi,jpj) ) |
|---|
| 393 | |
|---|
| 394 | ALLOCATE( ssmask(jpi,jpj), gru(jpi,jpj), grv(jpi,jpj), risfdep(jpi,jpj), & |
|---|
| 395 | hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), & |
|---|
| 396 | & uslpml(jpi,jpj), vslpml(jpi,jpj), wslpiml(jpi,jpj), wslpjml(jpi,jpj) ) |
|---|
| 397 | |
|---|
| 398 | ALLOCATE( uslp(jpi,jpj,jpk), vslp(jpi,jpj,jpk), wslpi(jpi,jpj,jpk), wslpj(jpi,jpj,jpk) ) |
|---|
| 399 | |
|---|
| 400 | ALLOCATE( akz(jpi,jpj,jpk), ah_wslp2(jpi,jpj,jpk), omlmask(jpi,jpj,jpk) ) |
|---|
| 401 | |
|---|
| 402 | ALLOCATE( ts_b(jpi,jpj,jpk,jpts), ts_a(jpi,jpj,jpk,jpts), rab_b(jpi,jpj,jpk,jpts) ) |
|---|
| 403 | |
|---|
| 404 | ALLOCATE( ua(jpi,jpj,jpk), va(jpi,jpj,jpk) ) |
|---|
| 405 | |
|---|
| 406 | IF(ln_tiled) ALLOCATE( zftu(jpi,jpj,jpk,jpts), zftv(jpi,jpj,jpk,jpts) ) ! arrays used for diagnostics ! tiled |
|---|
| 407 | |
|---|
| 408 | !-------------------------------------------------------------------------------------------- |
|---|
| 409 | ! 6. Calculate mbkt, mbku, mbkv, ice-shelf fields, e3t_n, e3w_n, gdept, r1_e1u, r1_e2v |
|---|
| 410 | |
|---|
| 411 | !$ACC KERNELS |
|---|
| 412 | mbkt(:,:) = 0 ; mbku(:,:) = 0 ; mbkv(:,:) = 0 |
|---|
| 413 | DO jk = 1, jpk |
|---|
| 414 | DO jj = 1, jpj |
|---|
| 415 | DO ji = 1, jpi |
|---|
| 416 | IF ( tmask(ji,jj,jk) == 1._wp ) mbkt(ji,jj) = jk |
|---|
| 417 | IF ( umask(ji,jj,jk) == 1._wp ) mbku(ji,jj) = jk |
|---|
| 418 | IF ( vmask(ji,jj,jk) == 1._wp ) mbkv(ji,jj) = jk |
|---|
| 419 | END DO |
|---|
| 420 | END DO |
|---|
| 421 | END DO |
|---|
| 422 | |
|---|
| 423 | DO jj = 1, jpj |
|---|
| 424 | DO ji = 1, jpi |
|---|
| 425 | mbkt(ji,jj) = MAX( mbkt(ji,jj) , 1) |
|---|
| 426 | mbku(ji,jj) = MAX( mbku(ji,jj) , 1) |
|---|
| 427 | mbkv(ji,jj) = MAX( mbkv(ji,jj) , 1) |
|---|
| 428 | END DO |
|---|
| 429 | END DO |
|---|
| 430 | |
|---|
| 431 | ! set ice-shelf fields for no ice ! |
|---|
| 432 | mikt(:,:) = 1 ; miku(:,:) = 1 ; mikv(:,:) = 1 ; risfdep(:,:) = 0._wp |
|---|
| 433 | gtui(:,:,:) = 0._wp ; gtvi(:,:,:) = 0._wp |
|---|
| 434 | |
|---|
| 435 | ! The following code is taken from dommsk: dom_msk |
|---|
| 436 | wmask (:,:,1) = tmask(:,:,1) |
|---|
| 437 | DO jk = 2, jpk ! interior values |
|---|
| 438 | wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1) |
|---|
| 439 | END DO |
|---|
| 440 | |
|---|
| 441 | ! The following code is taken from domvvl: dom_vvl_interpol for case 'W' interpolation |
|---|
| 442 | |
|---|
| 443 | IF ( ln_linssh ) THEN |
|---|
| 444 | e3w_n(:,:,:) = e3w_0(:,:,:) |
|---|
| 445 | |
|---|
| 446 | ELSE |
|---|
| 447 | e3w_n(:,:,1) = e3w_0(:,:,1) + e3t_n(:,:,1) - e3t_0(:,:,1) |
|---|
| 448 | ! - ML - The use of mask in this formula enables the special treatment of the last w-point without indirect adressing |
|---|
| 449 | DO jk = 2, jpk |
|---|
| 450 | e3w_n(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( e3t_n(:,:,jk-1) - e3t_0(:,:,jk-1) ) & |
|---|
| 451 | & + 0.5_wp * tmask(:,:,jk) * ( e3t_n(:,:,jk ) - e3t_0(:,:,jk ) ) |
|---|
| 452 | END DO |
|---|
| 453 | |
|---|
| 454 | END IF |
|---|
| 455 | |
|---|
| 456 | ! The following code is taken from domhgr: SUBROUTINE dom_hgr |
|---|
| 457 | e1e2t (:,:) = e1t(:,:) * e2t(:,:) ; r1_e1e2t(:,:) = 1._wp / e1e2t(:,:) |
|---|
| 458 | e1e2u (:,:) = e1u(:,:) * e2u(:,:) ; r1_e1e2u(:,:) = 1._wp / e1e2u(:,:) |
|---|
| 459 | e1e2v (:,:) = e1v(:,:) * e2v(:,:) ; r1_e1e2v(:,:) = 1._wp / e1e2v(:,:) |
|---|
| 460 | |
|---|
| 461 | e2_e1u(:,:) = e2u(:,:) / e1u(:,:) ; e1_e2v(:,:) = e1v(:,:) / e2v(:,:) |
|---|
| 462 | |
|---|
| 463 | ! The following code is taken from domvvl |
|---|
| 464 | DO jk = 1, jpk |
|---|
| 465 | DO jj = 1, jpj-1 |
|---|
| 466 | DO ji = 1, jpi-1 |
|---|
| 467 | !* from T- to U-point : hor. surface weighted mean |
|---|
| 468 | e3u_n(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e1e2u(ji,jj) & |
|---|
| 469 | & * ( e1e2t(ji ,jj) * ( e3t_n(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & |
|---|
| 470 | & + e1e2t(ji+1,jj) * ( e3t_n(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) |
|---|
| 471 | e3u_n(ji,jj,jk) = e3u_n(ji,jj,jk) + e3u_0(ji,jj,jk) |
|---|
| 472 | |
|---|
| 473 | !* from T- to V-point : hor. surface weighted mean |
|---|
| 474 | e3v_n(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e1e2v(ji,jj) & |
|---|
| 475 | & * ( e1e2t(ji,jj ) * ( e3t_n(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & |
|---|
| 476 | & + e1e2t(ji,jj+1) * ( e3t_n(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) |
|---|
| 477 | e3v_n(ji,jj,jk) = e3v_n(ji,jj,jk) + e3v_0(ji,jj,jk) |
|---|
| 478 | END DO |
|---|
| 479 | END DO |
|---|
| 480 | END DO |
|---|
| 481 | |
|---|
| 482 | ! there needs to be an lbc_lnk call here for e3u_n and e3v_n |
|---|
| 483 | e3u_n(jpi,:,:) = e3u_n(2,:,:) ; e3v_n(jpi,:,:) = e3v_n(2,:,:) |
|---|
| 484 | e3u_n(:,jpj,:) = e3u_n(:,2,:) ; e3v_n(:,jpj,:) = e3v_n(:,2,:) |
|---|
| 485 | |
|---|
| 486 | ! The following code is taken from domvvl: dom_vvl_init |
|---|
| 487 | |
|---|
| 488 | gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) ! reference to the ocean surface (used for MLD and light penetration) |
|---|
| 489 | gdepw_n(:,:,1) = 0.0_wp |
|---|
| 490 | gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) |
|---|
| 491 | |
|---|
| 492 | DO jk = 2, jpk ! vertical sum |
|---|
| 493 | DO jj = 1,jpj |
|---|
| 494 | DO ji = 1,jpi |
|---|
| 495 | ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt |
|---|
| 496 | ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) |
|---|
| 497 | ! ! 0.5 where jk = mikt |
|---|
| 498 | zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) |
|---|
| 499 | gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) |
|---|
| 500 | gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) & |
|---|
| 501 | & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) |
|---|
| 502 | gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) |
|---|
| 503 | |
|---|
| 504 | END DO |
|---|
| 505 | END DO |
|---|
| 506 | END DO |
|---|
| 507 | |
|---|
| 508 | DO jj = 1, jpj |
|---|
| 509 | DO ji = 1, jpi |
|---|
| 510 | r1_e1u(ji,jj) = 1._wp / e1u(ji,jj) |
|---|
| 511 | r1_e2v(ji,jj) = 1._wp / e2v(ji,jj) |
|---|
| 512 | END DO |
|---|
| 513 | END DO |
|---|
| 514 | |
|---|
| 515 | !-------------------------------------------------------------------------------------------- |
|---|
| 516 | ! 7. Initialise fields for calculation of isopycnal diffusion |
|---|
| 517 | |
|---|
| 518 | uslp (:,:,:) = 0._wp ; uslpml (:,:) = 0._wp |
|---|
| 519 | vslp (:,:,:) = 0._wp ; vslpml (:,:) = 0._wp |
|---|
| 520 | wslpi(:,:,:) = 0._wp ; wslpiml(:,:) = 0._wp |
|---|
| 521 | wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp |
|---|
| 522 | |
|---|
| 523 | akz(:,:,:) = 0._wp ; ah_wslp2(:,:,:) = 0._wp ! just to set values in the halo regions |
|---|
| 524 | |
|---|
| 525 | ts_b(:,:,:,:) = ts_n(:,:,:,:) ! set set ts on before time-step to the same values as the now timestep |
|---|
| 526 | |
|---|
| 527 | ! simplified choice of vertical diffusivity based on zdfphy |
|---|
| 528 | avt (:,:, 1 ) = 0._wp ; avt (:,:,jpk) = 0._wp |
|---|
| 529 | avt(:,:,2:jpkm1) = rn_avt0 |
|---|
| 530 | |
|---|
| 531 | ! simplified choice of horizontal diffusivity based on ldftra: ldf_tra_init |
|---|
| 532 | ! assumes Laplacian rather than bi-Laplacian isopycnal diffusion (inn = 1) l CASE nn_aht_ijk_t = 20 |
|---|
| 533 | zUfac = r1_2 *rn_Ud |
|---|
| 534 | |
|---|
| 535 | ! from ldfc1d_c2d:ldf_c2d case TRA with knn = 1 |
|---|
| 536 | DO jj = 1, jpj |
|---|
| 537 | DO ji = 1, jpi |
|---|
| 538 | ahtu(ji,jj,:) = zUfac * MAX( e1u(ji,jj), e2u(ji,jj) ) |
|---|
| 539 | ahtv(ji,jj,:) = zUfac * MAX( e1v(ji,jj), e2v(ji,jj) ) |
|---|
| 540 | END DO |
|---|
| 541 | END DO |
|---|
| 542 | !$ACC END KERNELS |
|---|
| 543 | |
|---|
| 544 | ! taken from dommsk |
|---|
| 545 | !$ACC KERNELS |
|---|
| 546 | !ssmask (:,:) = MAXVAL( tmask(:,:,:), DIM=3 ) ! Causes a run time failure on GPU, work around below |
|---|
| 547 | DO jj = 1, jpj |
|---|
| 548 | DO ji = 1, jpi |
|---|
| 549 | ssmask(ji,jj)=MAXVAL(tmask(ji,jj,:), DIM=1) |
|---|
| 550 | END DO |
|---|
| 551 | END DO |
|---|
| 552 | !$ACC END KERNELS |
|---|
| 553 | |
|---|
| 554 | !-------------------------------------------------------------------------------------------- |
|---|
| 555 | ! 8. Calculate the pressure gradients |
|---|
| 556 | |
|---|
| 557 | CALL eos_init |
|---|
| 558 | |
|---|
| 559 | ! IF istep == 1 calculate pressure gradients without the bias pressure correction (bpc) ; if istep == 2 calculate them with the bpc |
|---|
| 560 | |
|---|
| 561 | WRITE(0,*)"Stepping ",maxstep,"steps" |
|---|
| 562 | IF ( ln_zco ) THEN |
|---|
| 563 | write(0,*) ' Calling hpg_zco' |
|---|
| 564 | ELSE IF ( ln_zps ) THEN |
|---|
| 565 | write(0,*) ' Calling hpg_zps' |
|---|
| 566 | ELSE IF ( ln_sco ) THEN |
|---|
| 567 | write(0,*) ' Calling hpg_sco' |
|---|
| 568 | END IF |
|---|
| 569 | IF(ln_tiled) THEN |
|---|
| 570 | WRITE(0,*)" Tiled isopycnal diffusion",ji_len_3D_tile, jj_len_3D_tile, jk_len_3D_tile |
|---|
| 571 | ELSE |
|---|
| 572 | WRITE(0,*)" Non-Tiled isopycnal diffusion" |
|---|
| 573 | ENDIF |
|---|
| 574 | WRITE(0,*) |
|---|
| 575 | |
|---|
| 576 | step_time = TIMER() |
|---|
| 577 | DO istep = 1, maxstep |
|---|
| 578 | |
|---|
| 579 | !$ACC KERNELS |
|---|
| 580 | IF ( istep == maxstep ) THEN |
|---|
| 581 | ts_n(:,:,:,:) = ts_n(:,:,:,:) + ts_pcbias(:,:,:,:) |
|---|
| 582 | END IF |
|---|
| 583 | |
|---|
| 584 | ua(:,:,:) = 0.0_wp ; va(:,:,:) = 0.0_wp |
|---|
| 585 | |
|---|
| 586 | !$ACC END KERNELS |
|---|
| 587 | |
|---|
| 588 | et = TIMER() |
|---|
| 589 | 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 |
|---|
| 590 | eos_time = eos_time + (TIMER() - et) |
|---|
| 591 | |
|---|
| 592 | IF( ln_zps ) THEN |
|---|
| 593 | CALL zps_hde( istep, jpts, ts_n, mbku, mbkv, e3w_n, gdept_n, tmask, umask, vmask, & ! Partial steps: before horizontal gradient |
|---|
| 594 | & gtsu, gtsv, rhd, gru , grv, & ! of t, s, rd at the last ocean level |
|---|
| 595 | & zti, zhi, zri, ztj, zhj, zrj ) |
|---|
| 596 | END IF |
|---|
| 597 | |
|---|
| 598 | IF ( ln_zco ) THEN ! z-coordinate |
|---|
| 599 | et = TIMER() |
|---|
| 600 | CALL hpg_zco( grav, r1_e1u, r1_e2v, e3w_n, rhd, ua, va, zhpi ,zhpj) |
|---|
| 601 | hpg_zco_time = hpg_zco_time + (TIMER() - et) |
|---|
| 602 | |
|---|
| 603 | ELSE IF ( ln_zps ) THEN ! z-coordinate plus partial steps (interpolation) |
|---|
| 604 | et = TIMER() |
|---|
| 605 | CALL hpg_zps( grav, r1_e1u, r1_e2v, mbku, mbkv, gru, grv, e3w_n, rhd, ua, va, zhpi ,zhpj) |
|---|
| 606 | hpg_zps_time = hpg_zps_time + (TIMER() - et) |
|---|
| 607 | |
|---|
| 608 | ELSE IF ( ln_sco ) THEN |
|---|
| 609 | et = TIMER() |
|---|
| 610 | CALL hpg_sco( grav, r1_e1u, r1_e2v, e3w_n, rhd, gde3w_n, ln_linssh, ua, va, zhpi ,zhpj ) |
|---|
| 611 | hpg_sco_time = hpg_sco_time + (TIMER() - et) |
|---|
| 612 | |
|---|
| 613 | END IF |
|---|
| 614 | |
|---|
| 615 | !$ACC KERNELS |
|---|
| 616 | IF (istep == 1) THEN |
|---|
| 617 | hpgi(:,:,:) = ua(:,:,:) |
|---|
| 618 | hpgj(:,:,:) = va(:,:,:) |
|---|
| 619 | ELSE |
|---|
| 620 | hbpcgi(:,:,:) = hpgi(:,:,:) - ua(:,:,:) |
|---|
| 621 | hbpcgj(:,:,:) = hpgj(:,:,:) - va(:,:,:) |
|---|
| 622 | END IF |
|---|
| 623 | |
|---|
| 624 | DO jk = 2, jpk ! vertical sum |
|---|
| 625 | DO jj = 1,jpj |
|---|
| 626 | DO ji = 1,jpi |
|---|
| 627 | IF ( umask(ji,jj,jk) == 0. ) ua(ji,jj,jk) = fill |
|---|
| 628 | IF ( vmask(ji,jj,jk) == 0. ) va(ji,jj,jk) = fill |
|---|
| 629 | END DO |
|---|
| 630 | END DO |
|---|
| 631 | END DO |
|---|
| 632 | !$ACC END KERNELS |
|---|
| 633 | |
|---|
| 634 | !$ACC KERNELS |
|---|
| 635 | sum_ua=SUM(ua) |
|---|
| 636 | sum_va=SUM(va) |
|---|
| 637 | !$ACC END KERNELS |
|---|
| 638 | |
|---|
| 639 | !-------------------------------------------------------------------------------------------- |
|---|
| 640 | ! 9. Calculate isopycnal diffusion |
|---|
| 641 | |
|---|
| 642 | et = TIMER() |
|---|
| 643 | |
|---|
| 644 | CALL eos_rab_3d( ts_b, gdept_n, tmask, rab_b ) |
|---|
| 645 | |
|---|
| 646 | !$ACC KERNELS |
|---|
| 647 | rn2b(:,:,1) = 0._wp ; rn2b(:,:,jpk) = 0._wp |
|---|
| 648 | !$ACC END KERNELS |
|---|
| 649 | |
|---|
| 650 | CALL bn2 ( ts_b, rab_b, gdepw_n, gdept_n, e3w_n, wmask, rn2b ) |
|---|
| 651 | |
|---|
| 652 | CALL zdf_mxl( rn2b, e3w_n, gdept_n, gdepw_n, avt, wmask, ssmask, mbkt, nlb10, nmln, hmld, hmlp, hmlpt ) |
|---|
| 653 | |
|---|
| 654 | CALL ldf_slp( ln_zps, ln_isfcav, rhd, rn2b, tmask, umask, vmask, wmask, gru, grv, & |
|---|
| 655 | & e1t, e2t, r1_e1u, r1_e2v, e3u_n, e3v_n, e3w_n, gdept_n, gdepw_n, mbku, mbkv, & |
|---|
| 656 | & mikt, miku, mikv, nmln, hmlp, hmlpt, risfdep, & |
|---|
| 657 | & omlmask, uslpml, vslpml, wslpiml, wslpjml, uslp, vslp, wslpi, wslpj ) |
|---|
| 658 | |
|---|
| 659 | IF(.NOT.ln_tiled) THEN |
|---|
| 660 | CALL tra_ldf_iso( ln_traldf_msc, ln_traldf_blp, ln_traldf_lap, ln_zps, ln_isfcav, jpts, kpass, r2dt, & |
|---|
| 661 | & e1t, e1u, e1v, e2t, e2u, e2v, e1e2t, e1_e2v, e2_e1u, r1_e1e2t, mbku, mbkv, miku, mikv, & |
|---|
| 662 | & umask, vmask, wmask, e3t_n, e3u_n, e3v_n, e3w_n, wslpi, wslpj, uslp, vslp, & |
|---|
| 663 | & ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts_b , ts_b, & |
|---|
| 664 | & akz, ah_wslp2, ts_a ) |
|---|
| 665 | ELSE |
|---|
| 666 | DO jsk = 1, jpkm1, jk_len_3D_tile |
|---|
| 667 | tile % jsk = jsk |
|---|
| 668 | tile % jsk_2 = MAX(jsk, 2) |
|---|
| 669 | tile % jskm1 = MAX(jsk-1, 1) |
|---|
| 670 | tile % jek = MIN(jsk + jk_len_3D_tile - 1, jpkm1) |
|---|
| 671 | tile % jekp1 = MIN(tile % jek + 1, jpkm1) |
|---|
| 672 | tile % jskh = tile % jsk - nn_hls |
|---|
| 673 | tile % jekh = tile % jek + nn_hls |
|---|
| 674 | |
|---|
| 675 | DO jsj = nn_hls + 1, jpj - nn_hls, jj_len_3D_tile |
|---|
| 676 | tile % jsj = jsj |
|---|
| 677 | tile % jsjm1 = jsj-1 |
|---|
| 678 | tile % jej = MIN(jsj + jj_len_3D_tile - 1, jpj - nn_hls) |
|---|
| 679 | tile % jejp1 = tile % jej+1 |
|---|
| 680 | tile % jsjh = tile % jsj - nn_hls |
|---|
| 681 | tile % jejh = tile % jej + nn_hls |
|---|
| 682 | DO jsi = nn_hls + 1, jpi - nn_hls, ji_len_3D_tile |
|---|
| 683 | tile % jsi = jsi |
|---|
| 684 | tile % jsim1 = jsi-1 |
|---|
| 685 | tile % jei = MIN(jsi + ji_len_3D_tile - 1, jpi - nn_hls) |
|---|
| 686 | tile % jeip1 = tile % jei+1 |
|---|
| 687 | tile % jsih = tile % jsi - nn_hls |
|---|
| 688 | tile % jeih = tile % jei + nn_hls |
|---|
| 689 | |
|---|
| 690 | CALL tra_ldf_iso_tile( tile, ln_traldf_msc, ln_traldf_blp, ln_traldf_lap, ln_zps, ln_isfcav, jpts, kpass, r2dt, & |
|---|
| 691 | & e1t, e1u, e1v, e2t, e2u, e2v, e1e2t, e1_e2v, e2_e1u, r1_e1e2t, mbku, mbkv, miku, mikv, & |
|---|
| 692 | & umask, vmask, wmask, e3t_n, e3u_n, e3v_n, e3w_n, wslpi, wslpj, uslp, vslp, & |
|---|
| 693 | & ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts_b , ts_b, ts_a, & |
|---|
| 694 | & l_hst, l_ptr, akz, ah_wslp2, zftu, zftv ) |
|---|
| 695 | |
|---|
| 696 | END DO ! jis |
|---|
| 697 | END DO ! jjs |
|---|
| 698 | END DO ! jks |
|---|
| 699 | ENDIF |
|---|
| 700 | diffusion_time = diffusion_time + (TIMER() - et) |
|---|
| 701 | |
|---|
| 702 | !$ACC KERNELS |
|---|
| 703 | sum_ts_a=SUM(ts_a(:,:,:,1)) |
|---|
| 704 | !$ACC END KERNELS |
|---|
| 705 | |
|---|
| 706 | WRITE(0,*) 'istep = ', istep |
|---|
| 707 | WRITE(0,*) 'UA SUM ',sum_ua |
|---|
| 708 | WRITE(0,*) 'VA SUM ',sum_va |
|---|
| 709 | WRITE(0,*) 'TSA SUM ',sum_ts_a |
|---|
| 710 | WRITE(0,*) 'Step time ', TIMER() - step_time |
|---|
| 711 | WRITE(0,*) |
|---|
| 712 | |
|---|
| 713 | END DO ! istep |
|---|
| 714 | |
|---|
| 715 | step_time = TIMER() - step_time |
|---|
| 716 | |
|---|
| 717 | !-------------------------------------------------------------------------------------------- |
|---|
| 718 | ! 10. Write out the pressure gradient fields |
|---|
| 719 | |
|---|
| 720 | write(0,*) 'writing fields: after pseudo time-steps; fields are written at tracer lats and lons for simplicity for now ' |
|---|
| 721 | |
|---|
| 722 | !iostat = NF90_PUT_VAR(ncid, nav_lon_varid, nav_lon) |
|---|
| 723 | !iostat = NF90_PUT_VAR(ncid, nav_lat_varid, nav_lat) |
|---|
| 724 | !iostat = NF90_PUT_VAR(ncid, nav_lev_varid, nav_lev) |
|---|
| 725 | !iostat = NF90_PUT_VAR(ncid, hpgi_varid, hpgi) |
|---|
| 726 | !iostat = NF90_PUT_VAR(ncid, hpgj_varid, hpgj) |
|---|
| 727 | !iostat = NF90_PUT_VAR(ncid, hbpcgi_varid, hbpcgi) |
|---|
| 728 | !iostat = NF90_PUT_VAR(ncid, hbpcgj_varid, hbpcgj) |
|---|
| 729 | IF(ln_tiled) THEN |
|---|
| 730 | ! !iostat = NF90_PUT_VAR(ncid, ta_varid, ts_a(:,:,:,1) ) |
|---|
| 731 | ! !iostat = NF90_PUT_VAR(ncid, sa_varid, ts_a(:,:,:,2) ) |
|---|
| 732 | ENDIF |
|---|
| 733 | !IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR putting variable: ', TRIM(nf90_strerror(iostat)) |
|---|
| 734 | ! |
|---|
| 735 | ! Check Validation based on GPU PGI result |
|---|
| 736 | IF ( ln_zco ) THEN ! z-coordinate |
|---|
| 737 | WRITE(0,*)"ZCO: Difference in UA",sum_ua+409653003427.6933_wp |
|---|
| 738 | WRITE(0,*)"ZCO: Difference in VA",sum_va+407945220514.0470_wp |
|---|
| 739 | ELSE IF ( ln_zps ) THEN ! z-coordinate plus partial steps (interpolation) |
|---|
| 740 | WRITE(0,*)"ZPS: Difference in UA",sum_ua+409653003427.7428_wp |
|---|
| 741 | WRITE(0,*)"ZPS: Difference in VA",sum_va+407945220514.0659_wp |
|---|
| 742 | ELSE IF ( ln_sco ) THEN |
|---|
| 743 | WRITE(0,*)"SCO: Difference in UA",sum_ua+409653003427.7501_wp |
|---|
| 744 | WRITE(0,*)"SCO: Difference in VA",sum_va+407945220514.1797_wp |
|---|
| 745 | END IF |
|---|
| 746 | |
|---|
| 747 | #ifndef MANAGE |
|---|
| 748 | !$ACC UPDATE HOST (hbpcgi, hbpcgj) |
|---|
| 749 | #endif |
|---|
| 750 | !-------------------------------------------------------------------------------------------- |
|---|
| 751 | ! 9. Close the output file and stop |
|---|
| 752 | ! |
|---|
| 753 | iostat = NF90_PUT_VAR(ncid, hbpcgi_varid, hbpcgi) |
|---|
| 754 | iostat = NF90_PUT_VAR(ncid, hbpcgj_varid, hbpcgj) |
|---|
| 755 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR putting variable: ', TRIM(nf90_strerror(iostat)) |
|---|
| 756 | iostat = NF90_CLOSE(ncid) |
|---|
| 757 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR closing file: ', TRIM(nf90_strerror(iostat)) |
|---|
| 758 | |
|---|
| 759 | IF ( ln_zco ) THEN ! z-coordinate |
|---|
| 760 | WRITE(0,*)'hpg_zco CPU time [s] is: ', hpg_zco_time |
|---|
| 761 | ELSE IF ( ln_zps ) THEN ! z-coordinate plus partial steps (interpolation) |
|---|
| 762 | WRITE(0,*)'hpg_zps CPU time [s] is: ', hpg_zps_time |
|---|
| 763 | ELSE IF ( ln_sco ) THEN |
|---|
| 764 | WRITE(0,*)'hpg_sco CPU time [s] is: ', hpg_sco_time |
|---|
| 765 | END IF |
|---|
| 766 | IF( ln_zps ) WRITE(0,*)'zps_hde CPU time [s] is: ', zps_hde_time |
|---|
| 767 | WRITE(0,*)'eos_time CPU time [s] is: ', eos_time |
|---|
| 768 | WRITE(0,*)'isopycnal diffusion CPU time [s] is: ', diffusion_time |
|---|
| 769 | WRITE(0,*)'Total step time CPU time [s] is: ', step_time |
|---|
| 770 | write(0,*) 'finishing' |
|---|
| 771 | |
|---|
| 772 | END PROGRAM main_bpc_hpg |
|---|
