I think the CVMix user guide will grow out of this document, but for now it just seems like an easy way to share with MOM and MPAS users what I've done in POP to make calls into CVMix.
POP has a module named vertical_mix that is used to call either constant mixing, Richardson mixing, or KPP mixing. These three methods have their own modules: vmix_const, vmix_rich, and vmix_kpp. (It looks like the tidal mixing is in tidal_mixing.F90, which is used by KPP, but I haven't gotten that far yet.) Routines from vertical_mix are called from the baroclinic module, which is the main driver for computing baroclinic velocities and tracer fields. All of the changes I describe below were made in vertical_mix.F90, and I suspect that MPAS and MOM will want to make changes in a similar module.
The first goal for CVMix was to change POP to call CVMix's vmix_background and vmix_convective modules rather than POP's vmix_const module. This was a pretty straight forward process, requiring three steps.
Constant mixing requires three packages from CVMix:
use vmix_kinds_and_types use vmix_background use vmix_convection |
The first package is needed to set up all the derived data types:
type (vmix_data_type), allocatable, dimension(:) :: Vmix_vars type (vmix_global_params_type) :: Vmix_params type (vmix_bkgnd_params_type) :: Vmix_bkgnd_params type (vmix_conv_params_type) :: Vmix_conv_params |
In POP, if you want to use constant and convective mixing, the vertical_mix module reads in various convection-related namelist variables, allocates memory for VDC and VVC (what POP calls the vertical diffusivity and viscosity coefficients, respectively) and then calls init_vmix_const() in the vmix_const module. That routine reads in POP's &vmix_const namelist and set the arrays VDC and VVC to the scalar constants const_vdc and const_vvc, respectively.
allocate (VDC(nx_block,ny_block,km,1,nblocks_clinic), & VVC(nx_block,ny_block,km, nblocks_clinic)) |
Note that nblocks_clinic is the number of threads, nx_block and ny_block are the number of longitude / latitude points in each thread, respectively, and km is the number of vertical levels. So there are nx_block*ny_block*nblocks_clinic columns, each with km levels.
Obviously CVMix will not read in a POP namelist, so I moved the input into the vertical_mix module. I also needed to initialize several CVMix variables:
ncol = nx_block*ny_block*nblocks_clinic allocate(Vmix_vars(ncol)) do bid=1,nblocks_clinic do inx=1,nx_block do jny=1,ny_block colid = (bid-1)*nx_block*ny_block+(jny-1)*nx_block+inx Vmix_vars(colid)%nlev = KMT(inx, jny, bid) end do end do end do Prandtl = 10.0_r8 call vmix_init(Vmix_vars, Vmix_params, ncol, km, Prandtl) |
Vmix_vars(:) is the main CVMix data type. It is an array over each column -- as mentioned before, in POP the number of columns is the product of three of the dimensions of VDC and VVC --, colid is a POP-specific renumbering scheme to use natural ordering to get the (nx,ny,nblock) index down to a single number. Besides allocating memory for Vmix_vars, this routine also assigns a value to each Vmix_vars(:)%nlev, which is the number of ocean levels that particular column has. (In POP, this information is stored in the array KMT.) The last step is to call vmix_init, which is currently in the vmix_kinds_and_types module (on Coffee Talk, Linda Richman would probably point out that "vmix_init is neither a kind nor a type, discuss"). This is where memory for each component of Vmix_vars(:) gets allocated and a couple global parameters (maximum number of levels, which is needed for allocating some memory in the Vmix_bkgnd_params_type, and the Prandtl number, which Steve has pointed out should really be tied to a mixing scheme rather than saved as a 'global' parameter).
Once Vmix_vars is set up, we are ready to set up the background and convective mixing parameters / initialize variable values:
call vmix_init_bkgnd(Vmix_bkgnd_params, const_vvc, const_vdc) if (convection_itype.eq.convect_type_diff) then call vmix_init_conv(Vmix_conv_params, convect_diff, convect_visc) end if do bid=1,nblocks_clinic do inx=1,nx_block do jny=1,ny_block VVC(inx,jny,:, bid) = Vmix_bkgnd_params%static_visc(1,1) VDC(inx,jny,:,1,bid) = Vmix_bkgnd_params%static_diff(1,1) end do end do end do |
vmix_init_bkgnd allocates memory for Vmix_bkgnd_params%static_visc and %static_diff; because POP is sending scalar constants, these will both be 1x1 arrays. If you send vmix_init_bkgnd() 2D arrays (and the number of columns), it will allocate more memory for those two components. Similarly, vmix_init_conv sets the value for the convective mixing in Vmix_conv_params. The second block then populates the POP mixing coefficients with the values set in the CVMix routine.
At each timestep, the call to update the mixing coefficients is straight-forward (note that in POP this is called inside a threaded region and bid is the block number or thread ID):
! Pack into vmix_kind_type do k=1,km kp1 = min(k+1,km) call state(k ,kp1,TMIX(:,:,k ,1), TMIX(:,:,k ,2), & this_block, RHOOUT=RHOK_MAT(:,:,k,bid) ) call state(kp1,kp1,TMIX(:,:,kp1,1), TMIX(:,:,kp1,2), & this_block, RHOOUT=RHOKP_MAT(:,:,k,bid) ) end do do inx=1,nx_block do jny=1,ny_block colid = (bid-1)*nx_block*ny_block+(jny-1)*nx_block+inx nlev = Vmix_vars(colid)%nlev Vmix_vars(colid)%dens(:) = RHOK_MAT (inx,jny,1:nlev,bid) Vmix_vars(colid)%dens_lwr(:) = RHOKP_MAT(inx,jny,1:nlev,bid) Vmix_vars(colid)%visc_iface(2:nlev+1) = VVC(inx,jny,1:nlev, bid) Vmix_vars(colid)%diff_iface(2:nlev+1,1) = VDC(inx,jny,1:nlev,1,bid) end do end do call vmix_coeffs_bkgnd(Vmix_vars, Vmix_bkgnd_params, & (bid-1)*nx_block*ny_block+1, & bid*nx_block*ny_block) if (convection_itype.eq.convect_type_diff) then call vmix_coeffs_conv(Vmix_vars, Vmix_conv_params, & (bid-1)*nx_block*ny_block+1, & bid*nx_block*ny_block) end if ! Unpack from vmix_kind_type do inx=1,nx_block do jny=1,ny_block colid = (bid-1)*nx_block*ny_block+(jny-1)*nx_block+inx nlev = Vmix_vars(colid)%nlev VVC(inx,jny,1:nlev, bid) = Vmix_vars(colid)%visc_iface(2:nlev+1) VDC(inx,jny,1:nlev,1,bid) = Vmix_vars(colid)%diff_iface(2:nlev+1,1) end do end do |
The two calls to state() are made to get the density at level k and k+1 (note that this takes two calls instead of just using RHOK_MAT(:,:,k,:) and RHOK_MAT(:,:,k+1,:) because the density at level k is calculated for a parcel that is adiabatically displaced down to level k+1). These densities are needed to determine where to use the convective mixing.
The calls to vmix_coeffs_bkgnd() and vmix_coeffs_conv() actually update the mixing coefficients in Vmix_vars(:), and then the data needs to be unpacked back into the POP variables.