From 5be40550b7d2c3fb9dc080fe266893695c7566c6 Mon Sep 17 00:00:00 2001 From: He Yanchun Date: Tue, 3 Mar 2026 16:29:51 +0100 Subject: [PATCH] Tracmass implementation for the isopycnic model MICOM add MICOM project folder with configurations for reading model fields and grid handle outcropping isopycnal layers (zero thickness) handle north periodic boundary contion (north-folding) shift u/v flux staggering of the C-grid from west-south to north-east use model output vertical flux (w_explicit); or tracmass computed wflx (may still have error for some particles) add water parcel depth as a tracer (defined as depth of layer interface; not accurate to be improved) minor fixes of Tracmass code --- .gitignore | 4 + projects/MICOM/Makefile.prj | 22 +++ projects/MICOM/kill_zones.F90 | 86 ++++++++++ projects/MICOM/namelist_MICOM.in | 220 ++++++++++++++++++++++++++ projects/MICOM/read_field.F90 | 262 +++++++++++++++++++++++++++++++ projects/MICOM/setup_grid.F90 | 43 +++++ src/mod_calendar.F90 | 2 +- src/mod_error.F90 | 22 +-- src/mod_init.F90 | 5 + src/mod_loop.F90 | 12 +- src/mod_pos_tstep.F90 | 111 ++++++++++--- src/mod_seed.F90 | 6 +- src/mod_swap.F90 | 6 +- src/mod_tracerf.F90 | 24 +-- src/mod_tracers.F90 | 6 +- src/mod_vars.F90 | 8 + src/mod_vertvel.F90 | 30 +++- src/mod_write.F90 | 1 + 18 files changed, 811 insertions(+), 59 deletions(-) create mode 100644 projects/MICOM/Makefile.prj create mode 100644 projects/MICOM/kill_zones.F90 create mode 100644 projects/MICOM/namelist_MICOM.in create mode 100644 projects/MICOM/read_field.F90 create mode 100644 projects/MICOM/setup_grid.F90 diff --git a/.gitignore b/.gitignore index 12117f8..1732cf5 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,8 @@ # Build and statics +*.csv +*.mod +*.o +_build docs/_build docs/_static docs/_templates diff --git a/projects/MICOM/Makefile.prj b/projects/MICOM/Makefile.prj new file mode 100644 index 0000000..ad64900 --- /dev/null +++ b/projects/MICOM/Makefile.prj @@ -0,0 +1,22 @@ +# Makefile for BLOM run +#======================================================================== + +# Vertical velocities +#------------------------------------------------------------------------ +#fl01 = -Dw_2dim # Turn off vertical velocities. +#fl01 = -Dw_3dim # Compute 3D vertical velocities. +fl01 = -Dw_explicit # Explicit 3D vertical velocities. + + +# Time integration scheme +#------------------------------------------------------------------------ +#fl02 = -Dtime_analytical + +# Isopycnic model +#------------------------------------------------------------------------ +fl03 = -Disopycnic_model + +#======================================================================== + +ORM_FLAGS= -D$(PROJECT) \ +$(fl01)$(fl02)$(fl03) diff --git a/projects/MICOM/kill_zones.F90 b/projects/MICOM/kill_zones.F90 new file mode 100644 index 0000000..25193ea --- /dev/null +++ b/projects/MICOM/kill_zones.F90 @@ -0,0 +1,86 @@ +SUBROUTINE kill_zones + !========================================================================== + ! + ! Purpose + ! ------- + ! + ! Defines the limits of the domain. + ! If a trajectory is outside the domain the subroutine will identify it + ! with a flag (nend) + ! + ! nend = 0 is reserved to the time limit + ! nend = 1 is reserved to trajectories reaching the surface + ! + ! ========================================================================== + + + USE mod_traj + USE mod_domain + USE mod_tracervars + + IMPLICIT NONE + + INTEGER :: nexit, itrac, numexit + + nend = -1 + + SELECT CASE(exitType) + + ! Exit domain defined by the boxes [ienw, iene]x[jens, jenn] + ! in the namelist + CASE(1) + DO nexit = 1, 10 + IF(ienw(nexit) <= x1 .AND. x1 <= iene(nexit) .AND. & + jens(nexit) <= y1 .AND. y1 <= jenn(nexit) ) THEN + nend = nexit +1 + END IF + END DO + + ! Exit domain defined by the value of tracers in the namelist tracere + CASE(2) + + numexit = MINLOC(tracerchoice, DIM=1,MASK=(tracerchoice==999)) - 1 + + DO nexit = 1, numexit + itrac = tracerchoice(nexit) + + IF ( maxormin(nexit)*tracervalue(itrac)>= maxormin(nexit)*tracere(nexit) ) THEN + nend = nexit +1 + END IF + + END DO + + ! Exit domain defined by the value of tracers in the namelist tracere and domain + CASE(3) + + ! First thermodynamic killing zone + numexit = MINLOC(tracerchoice, DIM=1,MASK=(tracerchoice==999)) - 1 + + DO nexit = 1, numexit + itrac = tracerchoice(nexit) + + IF ( maxormin(nexit)*tracervalue(itrac)>= maxormin(nexit)*tracere(nexit) ) THEN + nend = nexit +1 + END IF + + END DO + + ! Next the geographical killing zone + DO nexit = 1, 10 + IF(ienw(nexit) <= x1 .AND. x1 <= iene(nexit) .AND. & + jens(nexit) <= y1 .AND. y1 <= jenn(nexit) ) THEN + nend = nexit +1 + numexit + END IF + END DO + + ! If the exit domained is defined in a different way: + ! - Read from a file + ! - Defined by curve, no linear shape, ... + ! This must be hard coded below. + + CASE(4) + PRINT*, 'Hard coded limits of the domain' + + END SELECT + +END SUBROUTINE kill_zones diff --git a/projects/MICOM/namelist_MICOM.in b/projects/MICOM/namelist_MICOM.in new file mode 100644 index 0000000..8103eac --- /dev/null +++ b/projects/MICOM/namelist_MICOM.in @@ -0,0 +1,220 @@ +&INIT_GRID_DESCRIPTION + ! Direction of the grid (1 - E->W, S->N, B->U) + griddir = -1,1,1 + + ! Does the netcdf data start with a zero index? + zeroindx = .FALSE. + + l_onestep = .TRUE. + + ! Describe filenames + physDataDir = '/projects/NS9560K/noresm/cases/NHISTfrc2_f09_tn14_20191025/ocn/hist/' + physPrefixForm = 'NHISTfrc2_f09_tn14_20191025.micom.hm.', + dateformat = 'YYYY-MM' + + fileSuffix = '.nc', ! suffix of files, e.g. ".nc" + + ! Describe variable names + ueul_name = 'uflx', ! name of resolved zonal velocity + veul_name = 'vflx', ! ---- " ------- meridional + w_name = 'wflx', ! name of diagnosed vertical velocity + usgs_name = '', ! name of eddy induced zonal velocity + vsgs_name = '', ! ---- " ------- meridional +/ + +&INIT_GRID_SIZE + ! Size of global grid + imt = 360, ! number of global i points + jmt = 384, ! number of global j points + km = 53, ! number of global k points + + iperio = 1, ! zonal periodic + jperio = 2, ! meridional boundary conditions + + topoDataDir = './', ! directory of grid data + + hgridFile = 'grid.nc', ! file with dx, dy etc + !dy_name = 'pdy', ! name of dy at T point + dyu_name = 'udy', ! dy at U point + !dx_name = 'pdx', ! name of dx at T point + dxv_name = 'vdx', ! dx at V point + + dzt_name = 'dp' ! name of 3D dp at T point + !dzu_name = 'dp' ! --- " --- dp at U point + !dzv_name = 'dp' ! --- " --- dp at V point + +/ + +&INIT_GRID_SUBDOMAIN + l_subdom = .FALSE. ! Initialise subdomain + imindom = 50 ! Eastern and western wall for subdomain + imaxdom = 100 + jmindom = 50 ! Southern and northern wall for subdomain + jmaxdom = 100 +/ + +&INIT_GRID_TIME + ngcm_step = 1, ! time between two time levels + ! if ngcm_unit = 5, then ngcm_step = 3 means 3 months + ngcm_unit = 5, ! units of ngcm + + iter = 120, ! number of subcycles between time levels stopped +/ + +&INIT_START_DATE + startsec = 0, ! start time + startmin = 0, + starthour = 0, + startday = 1, ! start date + startmon = 1, + startyear = 2000, + noleap = .TRUE., ! don't use leap years? + mon30day = .true. +/ + +&INIT_RUN_TIME + ! Loop run + loopYears = .TRUE., ! make a time loop to run + ! e.g. 1000 years using only 100 years of data + loopStartYear = 2000, ! year to start loop from + loopEndYear = 2000, ! last year of loop + + ! Verbose + log_level = 1 ! precise 1->3->5 verbose + + ! Information of run time + intrun = 121 ! number of time steps for run +/ + +&INIT_WRITE_TRAJ + write_frec = 1 ! 1 = write at time intervals of gcm datasets (each ints) + ! 2 = write at each time iteration + ! 3 = write each spatial grid-crossing + ! 4 = write all the time steps + ! 5 = write only start and end positions + !write_form = 1, + write_form = 0, + outDataDir = './output/MICOM/', + outDataFile = 'output', + + timeformat = 1, ! Format of the time array + ! 0 - tt / 1 - ts / 2 - YYYY, MM, DD, HH + l_compress = .FALSE. ! Compress _run.csv files on exceeding 2GB filesize using gzip +/ + +&INIT_SEEDING + nff = 1, ! = 1 run forward trajectories + ! = -1 run backwards + isec = 1, ! = 1 start on zonal cell wall + ! = 2 start on meridional cell wall + ! = 3 start on vertical wall + idir = 1, ! = 1 start only when flux > 0 + ! = -1 only when flux < 0 + nqua = 2, ! number of trajectories can be set by + ! = 1 constant number of particles in all seeding cells + ! set by partQuant (particles / gridcell) + ! = 2 Each particle reflects mass transport at seeding. + ! set by partQuant (m3/s or kg/s per particle) + ! = 3 Each particle reflects air/water mass/volume at seeding. + ! set by partQuant in m3 or kg per particle + nsdtraj = 10, ! number of trajectories to allocate per seeding cell (used for allocation if nqua > 1) + partquant = 5e6, ! particles/gridcell or m3s-1/particle or m3/particle + loneparticle = 0, ! start only one trajectory for debugging + seedtype = 1, ! = 1 seed using ist,jst,kst (below) + ! = 2 seed using file with start positions + + ! when seedtype = 1 +! gulf stream + ist1 = 30, ! seed in box in range i = [ist1,ist2] + ist2 = 35, + jst1 = 266, ! j = [jst1,jst2] + jst2 = 266, + kst1 = 1, ! k = [kst1,kst2] + kst2 = 53, + ! when seedtype = 2 + seeddir = '', ! directory for seed files + seedfile = '', ! name of seed file + + + seedtime = 1, ! = 1 seed only for specific times + ! = 2 seed using file with start positions + + ! when seedtime = 1 + tst1 = 1, ! seed only between tst1 and tst2 + tst2 = 5, + + ! when seedtime = 2 + timefile = '', ! name of seed file +/ + +&INIT_TRACERS + l_tracers = .True. ! Should we evaluate tracers along the trajectories? + + ! Tracer name (description) + tracername = 'temp','saln','depth' + ! Tracer unit (description) + tracerunit = 'degC','PSU','m' + + ! Name of the variable in the netcdf (if it's read) + tracervarname = 'temp','saln','dz' + + ! Action (read or compute) + traceraction = 'read','read','read' + + ! Minimum and maximum value of the binning + tracermin = -3, 32, -11000 + tracermax = 33, 38, 0 + + ! Dimension of the tracer + tracerdimension = 3D, 3D, 3D +/ + +&INIT_TRACERS_SEEDING + ! minimum and maximum value for seeding + ! same order as in tracername + !tracer0min = 2, 34.9, + !tracer0max = 4, 35.0, +/ + +&INIT_KILLZONES + timax = 75000 ! time limit for traj [days] + l_nosurface = .TRUE. + + exitType = 1 + ! If exitType = 1: killing zone defined by a box + ! IMPORTANT -> The maximum number of boxes than can be defined are 10 + ! ienw,iene are western and eastern i indices of killzone + ! jens,jenn are southern and northern j indices of killzone + ! ienw = 0, 500, + ! iene = 100, 1000, + ! jens = 200, 450, + ! jenn = 300, 470, + ! sets up two killzones, the first i = [0,100], j = [200,300] + ! the other i = [500,1000], j = [450,470] + ienw = 10,10, ! west indices of killzones + iene = 150,150, ! east indices of killzones + jens = 365, 365, ! south indices of killzones + jenn = 366, 366, ! north indices of killzones + + ! If exitType = 2: killing zone defined by tracer value + ! Tracerchoice defines which tracer is used for the killing zone + ! it could be different + !tracerchoice = 1, 2 + !tracere = 25, 35 + !maxormin = 1, -1 + ! This set up defines a killing zone by the 25C isotherm (maximum) + ! and the 35psu isohaline (minimum) + +/ + +&INIT_POSTPROCESS + l_psi = .FALSE. ! Compute stream Functions + dirpsi = -1,-1,-1,-1 ! Direction of integration +/ + +&INIT_ACTIVE + l_diffusion = .FALSE. ! Apply Stochastic diffusion along particle trajectories + Ah = 0 ! Horizontal diffusivity coefficient + Av = 0 ! Vertical diffusivity coefficient + +/ diff --git a/projects/MICOM/read_field.F90 b/projects/MICOM/read_field.F90 new file mode 100644 index 0000000..5589057 --- /dev/null +++ b/projects/MICOM/read_field.F90 @@ -0,0 +1,262 @@ +SUBROUTINE read_field + + !========================================================================== + ! + ! Purpose + ! ------- + ! + ! Read test model output to advect trajectories. + ! Will be called by the loop each time step. + ! + ! Method + ! ------ + ! + ! Read velocities and optionally some tracers from netCDF files and + ! update velocity fields for TRACMASS. + ! + ! Updates the variables: + ! uflux and vflux + ! ========================================================================== + + + USE mod_precdef + USE mod_param + USE mod_vel + USE mod_time + USE mod_grid + USE mod_getfile + USE mod_tracervars + USE mod_tracers + USE mod_calendar + USE mod_swap + + USE netcdf + + IMPLICIT none + + INTEGER :: itrac + INTEGER :: k + + REAL(DP), ALLOCATABLE, DIMENSION(:,:,:) :: tmp3d,dz3d,tmptracer + CHARACTER (len=200) :: fieldFile, dateprefix + REAL, PARAMETER :: missing_value=9.9692100E+36 + + ! Reassign the time index of uflux and vflux, dzt, dzdt, hs, ... + CALL swap_time() + + ! Data files + dateprefix = ' ' + + ALLOCATE(tmp3d(imt,jmt,km),dz3d(imt,jmt,km),tmptracer(imt,jmt,km)) + tmp3d(:,:,:) = 0.d0 + dz3d(:,:,:) = 0.d0 + tmptracer(:,:,:) = 0.d0 + + ! Reading 3-time step variables + ! In this case: dzt + ! =========================================================================== + IF (ints == 0) THEN + + ! 1 - Past + IF (loopYears) THEN + IF (l_onestep) nctstep = 1 + + dateprefix = filledFileName(dateFormat, prevYear, prevMon, prevDay) + + fieldFile = TRIM(physDataDir)//TRIM(physPrefixForm)//TRIM(dateprefix)//TRIM(fileSuffix) + dzt(:,:,km:1:-1,1) = get3DfieldNC(fieldFile,dzt_name,[imindom,jmindom,1,nctstep],[imt,jmt,km,1],'st') + WHERE (dzt(:,:,:,1) .eq. missing_value) dzt(:,:,:,1) = 0. + dzt(:,:,:,1) = dzt(:,:,:,1)/grav + + CALL count_ke(dzt(:,:,:,1),ke(:,:,1)) + !CALL write_ke(ke(:,:,1)) + END IF + + ! 2 - Present + nctstep = currMon + IF (l_onestep) nctstep = 1 + + dateprefix = filledFileName(dateFormat, currYear, currMon, currDay) + + fieldFile = TRIM(physDataDir)//TRIM(physPrefixForm)//TRIM(dateprefix)//TRIM(fileSuffix) + dzt(:,:,km:1:-1,2) = get3DfieldNC(fieldFile,dzt_name,[imindom,jmindom,1,nctstep],[imt,jmt,km,1],'st') + WHERE (dzt(:,:,:,2) .eq. missing_value) dzt(:,:,:,2) = 0. + dzt(:,:,:,2) = dzt(:,:,:,2)/grav ! dp -> dp/g + + CALL count_ke(dzt(:,:,:,2),ke(:,:,2)) + !CALL write_ke(ke(:,:,2)) + + END IF + + ! 3 - Future + IF (ints dp/g + WHERE (dzt(:,:,:,3) .eq. missing_value) dzt(:,:,:,3) = 0. + + CALL count_ke(dzt(:,:,:,2),ke(:,:,2)) + + END IF + + ! dzdt calculation + IF (ints == 0 .AND. ( loopYears .EQV..FALSE.)) THEN + dzdt(1:imt,1:jmt,:,2) = (dzt(1:imt,1:jmt,:,3) - dzt(1:imt,1:jmt,:,2))/tseas + ELSE IF (ints == intrun-1 .AND. ( loopYears .EQV..FALSE.)) THEN + dzdt(1:imt,1:jmt,:,2) = (dzt(1:imt,1:jmt,:,2) - dzt(1:imt,1:jmt,:,1))/tseas + ELSE + dzdt(1:imt,1:jmt,:,2) = 0.5*(dzt(1:imt,1:jmt,:,3) - dzt(1:imt,1:jmt,:,1))/tseas + END IF + + ! Reading 2-time step variables + ! In this case: velocities and tracers + ! =========================================================================== + + nctstep = currMon + IF (l_onestep) nctstep = 1 + + dateprefix = filledFileName(dateFormat, currYear, currMon, currDay) + + fieldFile = TRIM(physDataDir)//TRIM(physPrefixForm)//TRIM(dateprefix)//TRIM(fileSuffix) + + uvel(1:imt,1:jmt,km:1:-1) = get3DfieldNC(fieldFile, ueul_name,[imindom,jmindom,1,nctstep],[imt,jmt,km,1],'st') + IF (usgs_name/='') THEN + tmp3d(1:imt,1:jmt,km:1:-1) = get3DfieldNC(fieldFile, usgs_name,[imindom,jmindom,1,nctstep],[imt,jmt,km,1],'st') + uvel(1:imt,1:jmt,1:km) = uvel(1:imt,1:jmt,1:km) + tmp3d(1:imt,1:jmt,1:km) + END IF + WHERE (dzt(:,:,:,2) <= dzteps ) uvel(:,:,:) = 0. + + vvel(1:imt,0:jmt,km:1:-1) = get3DfieldNC(fieldFile, veul_name,[imindom,jmindom,1,nctstep],[imt,jmt+1,km,1],'st') + IF (vsgs_name/='') THEN + tmp3d(1:imt,0:jmt,km:1:-1) = get3DfieldNC(fieldFile, vsgs_name,[imindom,jmindom,1,nctstep],[imt,jmt+1,km,1],'st') + vvel(1:imt,0:jmt,1:km) = vvel(1:imt,0:jmt,1:km) + tmp3d(1:imt,0:jmt,1:km) + END IF + ! set vvel=0 if dztimtdom .OR. ib<1 .OR. ib>imtdom .OR. & - ja<1 .OR. ja>jmtdom .OR. jb<1 .OR. jb>jmtdom .OR. & - y0<0 .OR. y0>jmtdom .OR. y1<0 .OR. y1>jmtdom .OR. & - x0<0 .OR. x0>imtdom .OR. x1<0 .OR. x1>imtdom .OR. & - ka<1 .OR. ka>km .OR. kb<1 .OR. kb>km .OR. & - z0<0 .OR. z0>km .OR. z1<0 .OR. z1>km ) THEN + IF (ia<1 .OR. ia>imtdom .OR. ib<1 .OR. ib>imtdom .OR. & + ja<1 .OR. ja>jmtdom .OR. jb<1 .OR. jb>jmtdom .OR. & + y0<0 .OR. y0>jmtdom .OR. y1<0 .OR. y1>jmtdom .OR. & + x0<0 .OR. x0>imtdom .OR. x1<0 .OR. x1>imtdom .OR. & + ka<1 .OR. ka>km .OR. kb<1 .OR. kb>km .OR. & + z0<0 .OR. z0>km .OR. z1<0 .OR. z1>km ) THEN nerror = nerror + 1 errCode = 3 @@ -92,7 +92,7 @@ SUBROUTINE errorCheck(teststr,errCode) CALL write_error(errCode) - END IF + END IF CASE ('landError') @@ -192,20 +192,20 @@ SUBROUTINE write_error(errCode) CASE(0) ! Include time - tt in seconds WRITE(53,"(I8,3(',',F13.5),2(',',F20.5),1(', ',A100))") ntrac, xw, yw, zw, subvol, tt, & - ADJUSTL(errorType(errCode)) + TRIM(ADJUSTL(errorType(errCode))) RETURN CASE(1) ! Include time - Fraction ts WRITE(53,"(I8,3(',',F13.5),1(',',F20.5),1(',',F13.5),1(', ',A100))") ntrac, xw, yw, zw, subvol, ts, & - ADJUSTL(errorType(errCode)) + TRIM(ADJUSTL(errorType(errCode))) RETURN CASE(2) ! Include time - YYYY MM DD HH MM SS CALL tt_calendar(tt) WRITE(53,"(I8,3(',',F13.5),1(',',F20.5),(',',I5),3(',',I3),1(', ',A100))") ntrac, xw, yw, zw, & - subvol, dateYear, dateMon, dateDay, dateHour, ADJUSTL(errorType(errCode)) + subvol, dateYear, dateMon, dateDay, dateHour, TRIM(ADJUSTL(errorType(errCode))) RETURN END SELECT diff --git a/src/mod_init.F90 b/src/mod_init.F90 index dfe76e3..6a10148 100644 --- a/src/mod_init.F90 +++ b/src/mod_init.F90 @@ -128,6 +128,11 @@ SUBROUTINE init_alloc() ALLOCATE ( kmt(imt,jmt)) kmt(:,:) = km +#if defined isopycnic_model + ALLOCATE (ke(imt,jmt,nst)) + ke = 0 +#endif + ALLOCATE( dzt(imt,jmt,km,nst+1), dzu(imt,jmt,km,nst), dzv(imt,jmt,km,nst)) dzt(:,:,:,:) = 0.; dzu(:,:,:,:) = 0.; dzv(:,:,:,:) = 0. diff --git a/src/mod_loop.F90 b/src/mod_loop.F90 index fb7a1ce..ad9aec3 100755 --- a/src/mod_loop.F90 +++ b/src/mod_loop.F90 @@ -100,6 +100,16 @@ SUBROUTINE loop niter = trajectories(ntrac)%niter ts = DBLE(trajectories(ntrac)%nts) +#ifdef isopycnic_model + ! why nsm, not nsp? + DO WHILE (dzt(ib,jb,kb,nsm) .LT. dzteps) + kb = kb +1 + z1 = DBLE(kb) + trajectories(ntrac)%kb = kb + trajectories(ntrac)%z1 = z1 + END DO +#endif + ! Error code errCode = 0 @@ -137,7 +147,6 @@ SUBROUTINE loop intrpg = DMOD(ts,1.d0) ! -> gets the fractional part intrpr = 1.d0-intrpg - IF(intrpg.LT.0.d0 .OR. intrpg.GT.1.d0) THEN PRINT *,'* intrpg = ',intrpg PRINT *,'* intrpr = ',intrpr @@ -201,7 +210,6 @@ SUBROUTINE loop #ifndef time_analytical CALL update_bounce(ia, iam, ja, ka, x0, y0, z0) #endif - CALL cross_time(1,ia,ja,ka,x0,dse,dsw) ! zonal CALL cross_time(2,ia,ja,ka,y0,dsn,dss) ! meridional CALL cross_time(3,ia,ja,ka,z0,dsu,dsd) ! vertical diff --git a/src/mod_pos_tstep.F90 b/src/mod_pos_tstep.F90 index 8a4e46e..4909117 100755 --- a/src/mod_pos_tstep.F90 +++ b/src/mod_pos_tstep.F90 @@ -27,7 +27,7 @@ MODULE mod_pos USE mod_stream USE mod_grid, only: undef, imt, jmt, km, nsm, nsp, iperio, jperio - USE mod_vel, only: uflux, vflux, wflux, uu, um, vv, vm + USE mod_vel, only: uflux, vflux, wflux, uu, um USE mod_time, only: dtreg, intrpr, intrpg IMPLICIT none @@ -76,7 +76,7 @@ SUBROUTINE cross_time(ijk,ia,ja,ka,r0,sp,sn) ELSEIF (ijk .EQ. 3) THEN ii=ka -#if defined w_explicit +#if defined w_explicit uu = (intrpg*wflux(ia,ja,ka ,nsp) + intrpr*wflux(ia,ja,ka ,nsm)) um = (intrpg*wflux(ia,ja,ka-1,nsp) + intrpr*wflux(ia,ja,ka-1,nsm)) #else @@ -217,6 +217,10 @@ SUBROUTINE update_traj(ia,iam,ja,ka,ib,jb,kb,x0,y0,z0,x1,y1,z1) REAL(DP), INTENT(INOUT) :: x0, y0, z0 ! current position (referenced to the box) REAL(DP), INTENT(OUT) :: x1, y1, z1 ! future position (referenced to the box) +#if defined isopycnic_model + INTEGER :: ne = 0 !number of empty cells +#endif + ! New position scrivi = .FALSE. @@ -229,7 +233,7 @@ SUBROUTINE update_traj(ia,iam,ja,ka,ib,jb,kb,x0,y0,z0,x1,y1,z1) ! Eastward grid-cell exit IF (ds==dse) THEN scrivi=.FALSE. - uu = (intrpg*uflux(ia,ja,ka,nsp) + intrpr*uflux(ia ,ja,ka,nsm)) + uu = (intrpg*uflux(ia,ja,ka,nsp) + intrpr*uflux(ia,ja,ka,nsm)) IF (uu .GT. 0.d0) THEN ib=ia+1 @@ -353,9 +357,28 @@ SUBROUTINE update_traj(ia,iam,ja,ka,ib,jb,kb,x0,y0,z0,x1,y1,z1) uu = (intrpg*wflux(ka,nsp) + intrpr*wflux(ka,nsm)) #endif - IF (uu .GT. 0.d0) kb=ka+1 - - z1=DBLE(ka) + IF (uu .GT. 0.d0) THEN + kb = ka+1 + +#if defined isopycnic_model + if (intrpg .LT. dtmin) then + ne = max(ke(ia,ja,1),ke(ia,ja,2)) + else + ne = ke(ia,ja,2) + end if + IF(ka + ne >= KM-2) THEN !! ka or kb ??? + ka = KM-2 + kb = KM-1 + z1=DBLE(kb) + ELSE + z1=DBLE(ka) + END IF +#else + z1=DBLE(ka) +#endif + ELSE + z1=DBLE(ka) + END IF IF (kb==KM+1) THEN ! prevent particles to cross the boundaries kb=KM @@ -388,7 +411,6 @@ SUBROUTINE update_traj(ia,iam,ja,ka,ib,jb,kb,x0,y0,z0,x1,y1,z1) ! Downward grid-cell exit ELSE IF (ds==dsd) THEN - scrivi=.FALSE. CALL vertvel(ia, iam, ja, ka) #if defined w_explicit @@ -397,9 +419,27 @@ SUBROUTINE update_traj(ia,iam,ja,ka,ib,jb,kb,x0,y0,z0,x1,y1,z1) uu = (intrpg*wflux(ka-1,nsp) + intrpr*wflux(ka-1,nsm)) #endif - IF (uu .LT. 0.d0) kb=ka-1 - - z1=DBLE(ka-1) + IF (uu .LT. 0.d0) THEN + kb = ka-1 +#if defined isopycnic_model + IF (intrpg .LT. dtmin) THEN + ne = max(ke(ia,ja,1),ke(ia,ja,2)) + ELSE + ne = ke(ia,ja,2) + END IF + IF(kb==KM-2) THEN !K.Doos + ka = ka-ne-1 + kb = ka + z1=DBLE(ka) + ELSE + z1=DBLE(kb) + END IF +#else + z1=DBLE(kb) +#endif + ELSE + z1=DBLE(kb) + END IF CALL calc_pos(1,ia,ja,ka,x0,x1,ds) ! zonal position CALL calc_pos(2,ia,ja,ka,y0,y1,ds) ! meridional position @@ -450,10 +490,15 @@ SUBROUTINE update_traj(ia,iam,ja,ka,ib,jb,kb,x0,y0,z0,x1,y1,z1) jb = JMTdom - jmindom + 1 ia=ib ; ja=jb x0 = x1; y0 = y1 + ELSE IF( y1 + jmindom -1 == DBLE(JMTdom) .AND. jperio == 2) THEN + x1 = DBLE(IMTdom+imindom) - x1 -1 + ib = IDINT(x1) + 1 + jb = JMTdom - jmindom + 1 + ia=ib ; ja=jb + x0 = x1; y0 = y1 END IF END IF - ! Make sure that trajectory is inside ib,jb,kb box IF (x1 /= DBLE(IDINT(x1))) ib = IDINT(x1)+1 IF (y1 /= DBLE(IDINT(y1))) jb = IDINT(y1)+1 @@ -512,6 +557,10 @@ SUBROUTINE update_bounce(ia, iam, ja, ka, x0, y0, z0) ! Fluxes REAL(DP) :: uu, um ! Fluxes +#if defined isopycnic_model + INTEGER :: ne = 0 !number of empty cells +#endif + ! Temporal storage of indexes tmpia = ia tmpiam = iam @@ -524,7 +573,7 @@ SUBROUTINE update_bounce(ia, iam, ja, ka, x0, y0, z0) uu = (intrpg*uflux(ia,ja,ka,nsp) + intrpr*uflux(ia,ja,ka,nsm)) IF (uu .GT. 0.d0) THEN - ! Redifine the indexes + ! Redefine the indexes tmpiam = ia tmpia = ia + 1 IF ( (tmpia .EQ. IMT + 1) .AND. (iperio .EQ. 1)) tmpia = 1 @@ -535,7 +584,7 @@ SUBROUTINE update_bounce(ia, iam, ja, ka, x0, y0, z0) um = (intrpg*uflux(iam,ja,ka,nsp) + intrpr*uflux(iam,ja,ka,nsm)) IF (um .LT. 0.d0) THEN - ! Redifine the indexes + ! Redefine the indexes tmpia = iam tmpiam = tmpia - 1 IF( (tmpiam .EQ. 0) .AND. (iperio .EQ. 1) ) tmpiam = IMT @@ -549,7 +598,7 @@ SUBROUTINE update_bounce(ia, iam, ja, ka, x0, y0, z0) uu = (intrpg*vflux(ia,ja,ka,nsp) + intrpr*vflux(ia,ja,ka,nsm)) IF (uu .GT. 0.d0) THEN - ! Redifine the indexes + ! Redefine the indexes tmpja = ja + 1 END IF @@ -558,7 +607,7 @@ SUBROUTINE update_bounce(ia, iam, ja, ka, x0, y0, z0) um = (intrpg*vflux(ia,ja-1,ka,nsp) + intrpr*vflux(ia,ja-1,ka,nsm)) IF (um .LT. 0.d0) THEN - ! Redifine the indexes + ! Redefine the indexes tmpja = ja - 1 END IF @@ -575,8 +624,21 @@ SUBROUTINE update_bounce(ia, iam, ja, ka, x0, y0, z0) uu = (intrpg*wflux(ka ,nsp) + intrpr*wflux(ka ,nsm)) #endif IF (uu .GT. 0.d0) THEN - ! Redifine the indexes + ! Redefine the indexes +#if defined isopycnic_model + IF (intrpg .LT. dtmin) THEN + ne = max(ke(ia,ja,1),ke(ia,ja,2)) + ELSE + ne = ke(ia,ja,2) + END IF + IF(ka + ne >= KM-2) THEN + tmpka = KM-2 + ELSE tmpka = ka + 1 + END IF +#else + tmpka = ka + 1 +#endif END IF ELSE IF (z0 == DBLE(ka-1)) THEN @@ -589,9 +651,22 @@ SUBROUTINE update_bounce(ia, iam, ja, ka, x0, y0, z0) #endif IF (um .LT. 0.d0) THEN - ! Redifine the indexes + ! Redefine the indexes +#if defined isopycnic_model + IF (intrpg .LT. dtmin) THEN + ne = max(ke(ia,ja,1),ke(ia,ja,2)) + ELSE + ne = ke(ia,ja,2) + END IF + IF(ka-1==KM-2) THEN !K.Doos + tmpka = ka-ne-1 + ELSE tmpka = ka - 1 - END IF + END IF +#else + tmpka = ka - 1 +#endif + END IF END IF diff --git a/src/mod_seed.F90 b/src/mod_seed.F90 index 0c4d9ba..b912846 100644 --- a/src/mod_seed.F90 +++ b/src/mod_seed.F90 @@ -436,7 +436,7 @@ SUBROUTINE seed() PRINT *, ' vflux :', vflux (ib, jb, kb,nsm) PRINT *, ' subvol :', subvol STOP - ENDIF + END IF ! Determine start position for each particle ! -------------------------------------------------- @@ -588,7 +588,7 @@ SUBROUTINE seed() END IF CALL write_data('ini') - CALL write_data('run') + !CALL write_data('run') END DO kkkLoop END DO ijjLoop @@ -711,7 +711,7 @@ SUBROUTINE split_grid() ikt = jsg END IF - print*, num, ijt, ikt + print*, 'num, ijt, ikt:', num, ijt, ikt END IF diff --git a/src/mod_swap.F90 b/src/mod_swap.F90 index a3e1512..3a14da8 100644 --- a/src/mod_swap.F90 +++ b/src/mod_swap.F90 @@ -14,7 +14,7 @@ MODULE mod_swap !!------------------------------------------------------------------------------ USE mod_vel, only : uflux, vflux, wflux - USE mod_grid, only : dzt, dzu, dzv, dzdt, zstot, hs + USE mod_grid, only : dzt, dzu, dzv, dzdt, zstot, hs, ke USE mod_tracervars, only : l_tracers, l_swtraj, tracers, numtracers, tracertraj USE mod_time, only : nff @@ -51,7 +51,9 @@ SUBROUTINE swap_time() dzdt(:,:,:,1) = dzdt(:,:,:,2) - +#if defined isopycnic_model + ke(:,:,1) = ke(:,:,2) +#endif ! Scale factors and surface variables zstot(:,:,-1) = zstot(:,:,0) zstot(:,:, 0) = zstot(:,:,1) diff --git a/src/mod_tracerf.F90 b/src/mod_tracerf.F90 index e9b0b55..a728bf6 100644 --- a/src/mod_tracerf.F90 +++ b/src/mod_tracerf.F90 @@ -9,6 +9,8 @@ MODULE mod_tracerf !! !!------------------------------------------------------------------------------ + USE mod_precdef, only: DP + IMPLICIT NONE CONTAINS @@ -27,14 +29,14 @@ FUNCTION thermo_dens0(T,S) IMPLICIT NONE - REAL, INTENT(IN) :: T(:,:,:) ! Potential T [degC] - REAL, INTENT(IN) :: S(:,:,:) ! Practical S [PSU] + REAL(DP), INTENT(IN) :: T(:,:,:) ! Potential T [degC] + REAL(DP), INTENT(IN) :: S(:,:,:) ! Practical S [PSU] - REAL, ALLOCATABLE, DIMENSION (:,:,:) :: T68 + REAL(DP), ALLOCATABLE, DIMENSION (:,:,:) :: T68 - REAL, ALLOCATABLE, DIMENSION (:,:,:) :: thermo_dens0 + REAL(DP), ALLOCATABLE, DIMENSION (:,:,:) :: thermo_dens0 - REAL, ALLOCATABLE, DIMENSION (:,:,:) :: dens_temp + REAL(DP), ALLOCATABLE, DIMENSION (:,:,:) :: dens_temp INTEGER :: nx, ny, nz @@ -86,15 +88,15 @@ FUNCTION thermo_pt2ct(T,S) IMPLICIT NONE - REAL, INTENT(IN) :: T(:,:,:) ! Potential T [degC] - REAL, INTENT(IN) :: S(:,:,:) ! Absolut S [g/kg] + REAL(DP), INTENT(IN) :: T(:,:,:) ! Potential T [degC] + REAL(DP), INTENT(IN) :: S(:,:,:) ! Absolut S [g/kg] - REAL, ALLOCATABLE, DIMENSION (:,:,:) :: thermo_pt2ct + REAL(DP), ALLOCATABLE, DIMENSION (:,:,:) :: thermo_pt2ct - REAL, ALLOCATABLE, DIMENSION (:,:,:) :: ct_temp, xS, xS2, yT + REAL(DP), ALLOCATABLE, DIMENSION (:,:,:) :: ct_temp, xS, xS2, yT - REAL, PARAMETER :: gsw_sfac = 0.0248826675584615 - REAL, PARAMETER :: gsw_cp0 = 3991.86795711963 + REAL(DP), PARAMETER :: gsw_sfac = 0.0248826675584615 + REAL(DP), PARAMETER :: gsw_cp0 = 3991.86795711963 INTEGER :: nx, ny, nz diff --git a/src/mod_tracers.F90 b/src/mod_tracers.F90 index 19f4a13..4ba71d3 100644 --- a/src/mod_tracers.F90 +++ b/src/mod_tracers.F90 @@ -101,15 +101,15 @@ SUBROUTINE compute_tracer(tracname, var3d) ! Sigma 0 calculation (using T (degC) and S(g/kg)) IF (TRIM(tracname) == 'sigma0') THEN - var3d = REAL(thermo_dens0(REAL(tracers(1)%data(:,:,:,2),4), REAL(tracers(2)%data(:,:,:,2),4)),8) + var3d = thermo_dens0(tracers(1)%data(:,:,:,2), tracers(2)%data(:,:,:,2)) var3d = var3d - 1000.d0 ! Sigma 0 calculation (using T (K) and S(g/kg)) ELSE IF (TRIM(tracname) == 'sigma0_K') THEN - var3d = REAL(thermo_dens0(REAL(tracers(1)%data(:,:,:,2),4)-273.15, REAL(tracers(2)%data(:,:,:,2),4)),8) + var3d = thermo_dens0(tracers(1)%data(:,:,:,2)-273.15, tracers(2)%data(:,:,:,2)) var3d = var3d - 1000.d0 ! Conservative temperatue (CTo) calculation (using T (degC) and S(g/kg)) ELSE IF (TRIM(tracname) == 'CTo') THEN - var3d = REAL(thermo_pt2ct(REAL(tracers(1)%data(:,:,:,2),4), REAL(tracers(2)%data(:,:,:,2),4)),8) + var3d = thermo_pt2ct(tracers(1)%data(:,:,:,2), tracers(2)%data(:,:,:,2)) END IF END SUBROUTINE compute_tracer diff --git a/src/mod_vars.F90 b/src/mod_vars.F90 index 6fcc26a..6259853 100755 --- a/src/mod_vars.F90 +++ b/src/mod_vars.F90 @@ -67,6 +67,11 @@ MODULE mod_param REAL(DP), PARAMETER :: radian = pi/180.d0 REAL(DP), PARAMETER :: deg = radius*radian REAL(DP), PARAMETER :: tday = 24.d0 * 3600.d0 + +#ifdef isopycnic_model + REAL(DP), PARAMETER :: dzteps = 10.d0 ! rho*dz=10; dp ~ 100; rho~1e3 ->so dz~0.01m +#endif + ENDMODULE mod_param ! Seed Variables @@ -217,6 +222,9 @@ MODULE mod_grid INTEGER, ALLOCATABLE, DIMENSION(:,:) :: kmt, kmu, kmv +#if defined isopycnic_model + INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: ke ! total empty cells +#endif ! Sea level REAL(DP), ALLOCATABLE, DIMENSION(:,:,:) :: hs diff --git a/src/mod_vertvel.F90 b/src/mod_vertvel.F90 index f34da2e..52b1778 100755 --- a/src/mod_vertvel.F90 +++ b/src/mod_vertvel.F90 @@ -10,21 +10,25 @@ MODULE mod_vertvel !! !!------------------------------------------------------------------------------ - USE mod_vel, only : nsm, nsp, uflux, vflux, wflux, uu, um, vv, vm + USE mod_vel, only : nsm, nsp, uflux, vflux, wflux USE mod_time, only : intrpr, intrpg, tseas USE mod_grid +#ifdef isopycnic_model + USE mod_param, only : dzteps +#endif IMPLICIT NONE INTEGER :: k = 0 + REAL :: fnsm = 1.d0, fnsp =1.d0 + CONTAINS SUBROUTINE vertvel(ix, ixm, jy, kz) INTEGER :: ix, ixm, jy, kz - ! 1 #if defined w_2dim || w_explicit ! If 2D_w no w --// -- If 3D_w, w is read in the readfield @@ -33,17 +37,27 @@ SUBROUTINE vertvel(ix, ixm, jy, kz) #else kloop: DO k = 1, kz +#if defined isopycnic_model + IF ( dzt(ix,jy,k,nsm) <= dzteps ) fnsm = 0.d0 + IF ( dzt(ix,jy,k,nsp) <= dzteps ) fnsp = 0.d0 +#endif + IF (k> km - kmt(ix,jy)) THEN - wflux(k, 1) = wflux(k-1, 1) - & - ( uflux(ix,jy,k,1) - uflux(ixm,jy,k,1) + vflux(ix,jy,k,1) - vflux(ix,jy-1,k,1) ) & - - dzdt(ix,jy,k,1)*dxdy(ix,jy) + wflux(k,nsm) = wflux(k-1,nsm) - & + ( uflux(ix,jy,k,nsm) - uflux(ixm,jy,k,nsm) + vflux(ix,jy,k,nsm) - vflux(ix,jy-1,k,nsm) ) & + - dzdt(ix,jy,k,nsm)*dxdy(ix,jy)*fnsm - wflux(k, 2) = wflux(k-1, 2) - & - ( uflux(ix,jy,k,2) - uflux(ixm,jy,k,2) + vflux(ix,jy,k,2) - vflux(ix,jy-1,k,2) ) & - - dzdt(ix,jy,k,2)*dxdy(ix,jy) + wflux(k,nsp) = wflux(k-1,nsp) - & + ( uflux(ix,jy,k,nsp) - uflux(ixm,jy,k,nsp) + vflux(ix,jy,k,nsp) - vflux(ix,jy-1,k,nsp) ) & + - dzdt(ix,jy,k,nsp)*dxdy(ix,jy)*fnsp ELSE wflux(k,:) = 0.d0 END IF +#if defined isopycnic_model + ! reset + fnsm = 1.d0 + fnsp = 1.d0 +#endif END DO kloop #endif diff --git a/src/mod_write.F90 b/src/mod_write.F90 index 1ba8c5b..e89692c 100755 --- a/src/mod_write.F90 +++ b/src/mod_write.F90 @@ -222,6 +222,7 @@ SUBROUTINE write_data(sel) ( write_frec == 2 .AND. ABS(tss-DBLE(INT(tss)))<1e-11 .AND. ints == NINT(ts)) .OR. & ( write_frec == 3 .AND. (.not.scrivi .OR. boxface>0) ) .OR. & ( write_frec == 4 ) .OR. & + ( write_frec == 5 .AND. (trajectories(ntrac)%niter == niter-1) .OR. (.not.scrivi .OR. boxface>0) ) .OR. & ( write_frec == 2 .AND. tt == 0.d0)) THEN ! If postprocessing is activated