Description |
In this update, several files were modified, so the PIO library works in applications with tidal and river forcing.
- The reading routines nf_fread2d.F, nf_fread3d.F, nf_fread4d.F, nf_fread2d_bry.F, and nf_fread3d_bry.F were updated to compute the checksum values in a compact way. For example, in nf_fread3d.F, we have:
IF (Lchecksum) THEN
Npts=(Imax-Imin+1)*(Jmax-Jmin+1)*(UBk-LBk+1)
IF (.not.associated(Cwrk)) allocate ( Cwrk(Npts) )
Cwrk=PACK(Adat(Imin:Imax, Jmin:Jmax, LBk:UBk), .TRUE.)
CALL get_hash (Cwrk, Npts, checksum, .TRUE.)
IF (associated(Cwrk)) deallocate (Cwrk)
END IF
- Corrected routine pio_nfread2d in nf_fread2d.F module to work processing fields that need regridding.
- Restructured get_ngfld.F and get_ngfldr.F to clearly distinguish when processing input arrays of rank 1, 2, and 3. The PIO library expects the start and total arguments to NetCDF reading with congruent ranks. In addition, those routines have a new logical variable argument, recordless, to distinguish when processing variables with time records or not:
!
! Tidal Period.
!
IF (LprocessTides(ng)) THEN
IF (iic(ng).eq.0) THEN
CALL get_ngfld (ng, iNLM, idTper, TIDE(ng)%ncid, &
# if defined PIO_LIB && defined DISTRIBUTE
& TIDE(ng)%pioFile, &
# endif
& 1, TIDE(ng), recordless, update(1), &
& 1, MTC, 1, 1, 1, NTC(ng), 1, &
& TIDES(ng) % Tperiod)
IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN
END IF
END IF
#endif
All the calls to get_ngfld in get_data.F, rp_get_data.F, and tl_get_data.F were updated. Similarly, all the calls to get_ngfldr in ad_get_data.F.
Many thanks to Dave Robertson for bringing this issue to my attention.
|
Description |
This update enhances several routines:
- The routine ad_main3d is rewritten to include the nesting logic. The 4D-Var nesting scheme are working and will be released sometime next year. The complicated nested 4D-Var algorithms were developed in collaboration with Andy Moore. It was a herculean and tricky effort involving lots of testing.
- For consistency, the time-stepping counter is updated at the bottom of main3d, tl_main3d, and ad_main3d. For example, in main3d.F we have:
!
!-----------------------------------------------------------------------
! Advance time index and time clock.
!-----------------------------------------------------------------------
!
DO ig=1,GridsInLayer(nl)
ng=GridNumber(ig,nl)
iic(ng)=iic(ng)+1
time(ng)=time(ng)+dt(ng)
step_counter(ng)=step_counter(ng)-1
CALL time_string (time(ng), time_code(ng))
END DO
It cleaned the awkward logic at the bottom of initial, tl_initial, and ad_initial. This change was required for a more straightforward nested 4D-Var design.
- Adds AD_OUTPUT_STATE option to process and write out the whole adjoint solution. Due to the exact discrete adjoint, the predictor/corrector time-stepping scheme with multiple time levels, pieces of the adjoint solution are in two-time levels and need to be added in the ad_*_sol arrays for output purposes. We need this capability in the stochastic algorithms and the ROMS-JEDI interface.
Thus, changes are made to ad_ini_fields, ad_main2d, ad_main3d, and ad_step3d_uv. The two pieces for ad_u_sol and ad_v_sol are computed in ad_step3d_uv (index nnew). The other piece is computed in ad_pre_step3d, and added in ad_out_fields:
!
! Add a second piece of the 3D adjoint momentum solution into IO arrays.
! The first is loaded in "ad_step3d_uv".
!
DO j=JstrB,JendB
IF (j.ge.JstrM) THEN
DO k=1,N(ng)
DO i=IstrB,IendB
ad_v_sol(i,j,k)=ad_v_sol(i,j,k)+ad_v(i,j,k,nnew)
END DO
END DO
END IF
DO k=1,N(ng)
DO i=IstrM,IendB
ad_u_sol(i,j,k)=ad_u_sol(i,j,k)+ad_u(i,j,k,nnew)
END DO
END DO
END DO
- The time index logic in ad_wrt_his was slightly modified to account for the change of location of iic(ng) counter update in ad_main3d:
!
! Determine time index to write. The "nout" index is updated to the
! version of "ad_main3d" that updates the "iic" counter at the bottom.
! Therefore, we need to change the conditional "iic(ng).ne.ntend(ng)"
! to "iic(ng).gt.ntend(ng)" to get identical solutions.
!
# ifdef SOLVE3D
kout=kstp(ng)
# else
kout=kstp(ng)
# endif
# if defined WEAK_CONSTRAINT
kfout=2
# endif
# ifdef SOLVE3D
IF (iic(ng).gt.ntend(ng)) THEN
nout=nnew(ng)
# ifdef AD_OUTPUT_STATE
LwrtState3d(ng)=.FALSE.
# endif
ELSE
# ifdef AD_OUTPUT_STATE
LwrtState3d(ng)=.TRUE.
# endif
nout=nstp(ng)
END IF
# endif
- The jedi_roms.h driver was modified to facilitate the LinearModel Class for the ROMS-JEDI interface. OOPS will advance the ROMS kernels by smaller intervals, usually a single timestep. The strategy here is different from that used for coupling because ROMS delayed output. The delayed output last half timestep will affect the OOPS trajectory logic needed to save the NLM background fields required to linearize the TLM and ADM kernels. Therefore, we need to address some technical issues for the correct TL and AD state vector to exchange between ROMS and OOPS.
- Enhanced function netcdf_inq_var in mod_netcdf.F to include more information about file global dimensions: names, ID, and size.
- Updated the lateral boundary conditions for the vertical mixing coefficients in the turbulent parameterizations schemes. For example, in gls_corstep.F, we have:
DO k=0,N(ng)
IF (DOMAIN(ng)%Western_Edge(tile)) THEN
DO j=Jstr,Jend
DO itrc=1,NAT
Akt(Istr-1,j,k,itrc)=Akt(Istr,j,k,itrc)
END DO
Akv(Istr-1,j,k)=Akv(Istr,j,k)
END DO
END IF
IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
DO j=Jstr,Jend
DO itrc=1,NAT
Akt(Iend-1,j,k,itrc)=Akt(Iend,j,k,itrc)
END DO
Akv(Iend-1,j,k)=Akv(Iend,j,k)
END DO
END IF
IF (DOMAIN(ng)%Southern_Edge(tile)) THEN
DO i=Istr,Iend
DO itrc=1,NAT
Akt(i,Jstr-1,k,itrc)=Akt(i,Jstr,k,itrc)
END DO
Akv(i,Jstr-1,k)=Akv(i,Jstr,k)
END DO
END IF
IF (DOMAIN(ng)%Northern_Edge(tile)) THEN
DO i=Istr,Iend
DO itrc=1,NAT
Akt(i,Jend+1,k,itrc)=Akt(i,Jend,k,itrc)
END DO
Akv(i,Jend+1,k)=Akv(i,Jend,k)
END DO
END IF
IF (DOMAIN(ng)%SouthWest_Corner(tile)) THEN
DO itrc=1,NAT
Akt(Istr-1,Jstr-1,k,itrc)=0.5_r8* &
& (Akt(Istr ,Jstr-1,k,itrc)+ &
& Akt(Istr-1,Jstr ,k,itrc))
END DO
Akv(Istr-1,Jstr-1,k)=0.5_r8* &
& (Akv(Istr ,Jstr-1,k)+ &
& Akv(Istr-1,Jstr ,k))
END IF
IF (DOMAIN(ng)%SouthEast_Corner(tile)) THEN
DO itrc=1,NAT
Akt(Iend+1,Jstr-1,k,itrc)=0.5_r8* &
& (Akt(Iend ,Jstr-1,k,itrc)+ &
& Akt(Iend+1,Jstr ,k,itrc))
END DO
Akv(Iend+1,Jstr-1,k)=0.5_r8* &
& (Akv(Iend ,Jstr-1,k)+ &
& Akv(Iend+1,Jstr ,k))
END IF
IF (DOMAIN(ng)%NorthWest_Corner(tile)) THEN
DO itrc=1,NAT
Akt(Istr-1,Jend+1,k,itrc)=0.5_r8* &
& (Akt(Istr ,Jend+1,k,itrc)+ &
& Akt(Istr-1,Jend ,k,itrc))
END DO
Akv(Istr-1,Jend+1,k)=0.5_r8* &
& (Akv(Istr ,Jend+1,k)+ &
& Akv(Istr-1,Jend ,k))
END IF
IF (DOMAIN(ng)%NorthEast_Corner(tile)) THEN
DO itrc=1,NAT
Akt(Iend+1,Jend+1,k,itrc)=0.5_r8* &
& (Akt(Iend ,Jend+1,k,itrc)+ &
& Akt(Iend+1,Jend ,k,itrc))
END DO
Akv(Iend+1,Jend+1,k)=0.5_r8* &
& (Akv(Iend ,Jend+1,k)+ &
& Akv(Iend+1,Jend ,k))
END IF
END DO
The parameterization schemes compute the vertical mixing at the interior points, and we have boundary conditions for only periodic applications. This issue was discovered in the ROMS-JEDI interface. This change introduces minimal roundoff changes to the nonlinear model solution, which affects the TLM, RPM, ADM solutions.
- Corrected bug in def_avg_pio that gives error when synchronizing a PIO file to disk. Many thanks to Dave Robertson for bringing this issue to my attention.
|