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