MODULE dynspg_fsc !!====================================================================== !! *** MODULE dynspg_fsc *** !! Ocean dynamics: surface pressure gradient trend !!====================================================================== #if ( defined key_dynspg_fsc && ! defined key_autotasking ) || defined key_esopa !!---------------------------------------------------------------------- !! 'key_dynspg_fsc' free surface cst volume !! NOT 'key_autotasking' k-j-i loop (vector opt.) !!---------------------------------------------------------------------- !! dyn_spg_fsc : update the momentum trend with the surface pressure !! gradient in the free surface constant volume case !! with vector optimization !!---------------------------------------------------------------------- !! * Modules used USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE trdtra_oce ! ocean active tracer trend USE trddyn_oce ! ocean active dynamics trend USE zdf_oce ! ocean vertical physics USE in_out_manager ! I/O manager USE phycst ! physical constants USE ocesbc ! ocean surface boundary condition USE flxrnf ! ??? USE sol_oce ! solver variables USE solpcg ! preconditionned conjugate gradient solver USE solsor ! Successive Over-relaxation solver USE solfet ! FETI solver USE obc_oce ! Lateral open boundary condition USE obcdyn ! ocean open boundary condition (obc_dyn routines) USE obcvol ! ocean open boundary condition (obc_vol routines) USE lib_mpp ! ??? USE lbclnk ! ??? USE cla_dynspg ! ??? IMPLICIT NONE PRIVATE !! * Accessibility PUBLIC dyn_spg_fsc ! routine called by step.F90 !! * Shared module variables LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_fsc = .TRUE. !: free surface constant volume flag !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! OPA 9.0 , LODYC-IPSL (2003) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE dyn_spg_fsc( kt, kindic ) !!---------------------------------------------------------------------- !! *** routine dyn_spg_fsc *** !! !! ** Purpose : Compute the now trend due to the surface pressure !! gradient in case of free surface formulation with a constant !! ocean volume add it to the general trend of momentum equation. !! !! ** Method : Free surface formulation. The surface pressure gradient !! is given by: !! spgu = 1/rau0 d/dx(ps) = 1/e1u di( etn + rnu btda ) !! spgv = 1/rau0 d/dy(ps) = 1/e2v dj( etn + rnu btda ) !! where etn is the free surface elevation and btda is the after !! of the free surface elevation !! -1- compute the after sea surface elevation from the cinematic !! surface boundary condition: !! zssha = sshb + 2 rdt ( wn - emp ) !! Time filter applied on now sea surface elevation to avoid !! the divergence of two consecutive time-steps and swap of free !! surface arrays to start the next time step: !! sshb = sshn + atfp * [ sshb + zssha - 2 sshn ] !! sshn = zssha !! -2- evaluate the surface presure trend (including the addi- !! tional force) in three steps: !! a- compute the right hand side of the elliptic equation: !! gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ] !! where (spgu,spgv) are given by: !! spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ] !! - grav 2 rdt hu /e1u di[sshn + emp] !! spgv = vertical sum[ e3v (vb+ 2 rdt va) ] !! - grav 2 rdt hv /e2v dj[sshn + emp] !! and define the first guess from previous computation : !! zbtd = btda !! btda = 2 zbtd - btdb !! btdb = zbtd !! b- compute the relative accuracy to be reached by the !! iterative solver !! c- apply the solver by a call to sol... routine !! -3- compute and add the free surface pressure gradient inclu- !! ding the additional force used to stabilize the equation. !! !! ** Action : - Update (ua,va) with the surf. pressure gradient trend !! - Save the trends in (utrd,vtrd) ('key_diatrends') !! !! References : !! Roullet and Madec 1999, JGR. !! !! History : !! ! 98-05 (G. Roullet) Original code !! ! 98-10 (G. Madec, M. Imbard) release 8.2 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module !! ! 02-11 (C. Talandier, A-M Treguier) Open boundaries !!--------------------------------------------------------------------- !! * Arguments INTEGER, INTENT( in ) :: kt ! ocean time-step index INTEGER, INTENT( out ) :: kindic ! solver convergence flag ! if the solver doesn t converge ! the flag is < 0 !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & z2dt, z2dtg, zraur, znugdt, & ! temporary scalars znurau, zssha, zspgu, zspgv, & ! " " zegu, zegv, zgcb, zbtd, & ! " " ztdgu, ztdgv, zgwgt ! " " !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_spg_fsc : surface pressure gradient trend' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ (free surface constant volume case)' ! set to zero free surface specific arrays spgu(:,:) = 0.e0 ! surface pressur gradient (i-direction) spgv(:,:) = 0.e0 ! surface pressur gradient (j-direction) ENDIF ! 0. Local constant initialization ! -------------------------------- ! time step: leap-frog z2dt = 2. * rdt ! time step: Euler if restart from rest IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! coefficients z2dtg = grav * z2dt zraur = 1. / rauw znugdt = rnu * grav * z2dt znurau = znugdt * zraur IF( lk_mpp ) THEN ! Mpp : export boundary values of to neighboring processors !!bug : I don t understand why this only in mpp???? ==> Can be suppressed, no? CALL lbc_lnk( ua, 'U', -1. ) CALL lbc_lnk( va, 'V', -1. ) ENDIF ! 1. Surface pressure gradient (now) ! ---------------------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zspgu = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj) zspgv = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj) zegu = + znurau * ( emp (ji+1,jj) - emp (ji,jj) ) / e1u(ji,jj) zegv = + znurau * ( emp (ji,jj+1) - emp (ji,jj) ) / e2v(ji,jj) spgu(ji,jj) = zspgu + zegu spgv(ji,jj) = zspgv + zegv #if defined key_trddyn ! save the surface pressure gradient trends utrd2(ji,jj,1) = zspgu vtrd2(ji,jj,1) = zspgv ! save the first part of the additional force trend utrd2(ji,jj,2) = zegu vtrd2(ji,jj,2) = zegv #endif END DO END DO ! 2. Add the surface pressure trend to the general trend ! ------------------------------------------------------ DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) #if defined key_trddyn ! save the surface pressure gradient trend for diagnostics utrd(ji,jj,jk,8) = spgu(ji,jj) vtrd(ji,jj,jk,8) = spgv(ji,jj) #endif END DO END DO END DO ! 1. Evaluate the masked next velocity ! ------------------------------------ ! (effect of the additional force not included) DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk) va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk) END DO END DO END DO #if defined key_obc ! Update velocities on each open boundary with the radiation algorithm CALL obc_dyn( kt ) ! Correction of the barotropic componant velocity to control the volume of the system CALL obc_vol( kt ) #endif #if defined key_orca_r2 IF( n_cla == 1 ) CALL dyn_spg_cla( kt ) ! Cross Land Advection (update (ua,va)) #endif ! 2. compute the next vertically averaged velocity ! ------------------------------------------------ ! (effect of the additional force not included) ! initialize to zero DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. spgu(ji,jj) = 0.e0 spgv(ji,jj) = 0.e0 END DO END DO ! vertical sum !CDIR NOLOOPCHG IF( lk_vopt_loop ) THEN ! vector opt., forced unroll DO jk = 1, jpkm1 DO ji = 1, jpij spgu(ji,1) = spgu(ji,1) + fse3u(ji,1,jk) * ua(ji,1,jk) spgv(ji,1) = spgv(ji,1) + fse3v(ji,1,jk) * va(ji,1,jk) END DO END DO ELSE ! No vector opt. DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = 2, jpim1 spgu(ji,jj) = spgu(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk) spgv(ji,jj) = spgv(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk) END DO END DO END DO ENDIF ! transport: multiplied by the horizontal scale factor DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj) spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj) END DO END DO ! Boundary conditions on (spgu,spgv) CALL lbc_lnk( spgu, 'U', -1. ) CALL lbc_lnk( spgv, 'V', -1. ) ! 3. Right hand side of the elliptic equation and first guess ! ----------------------------------------------------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! Divergence of the after vertically averaged velocity zgcb = spgu(ji,jj) - spgu(ji-1,jj) & + spgv(ji,jj) - spgv(ji,jj-1) gcb(ji,jj) = gcdprc(ji,jj) * zgcb ! First guess of the after barotropic transport divergence zbtd = gcx(ji,jj) gcx (ji,jj) = 2. * zbtd - gcxb(ji,jj) gcxb(ji,jj) = zbtd END DO END DO ! 4. Relative precision (computation on one processor) ! --------------------- rnorme =0. DO jj = 1, jpj DO ji = 1, jpi zgwgt = gcdmat(ji,jj) * gcb(ji,jj) rnorme = rnorme + gcb(ji,jj) * zgwgt END DO END DO IF( lk_mpp ) CALL mpp_sum( rnorme ) ! sum over the global domain epsr = eps * eps * rnorme ncut = 0 ! if rnorme is 0, the solution is 0, the solver isn't called IF( rnorme == 0.e0 ) THEN gcx(:,:) = 0.e0 res = 0.e0 niter = 0 ncut = 999 ENDIF ! 5. Evaluate the next transport divergence ! ----------------------------------------- ! Iterarive solver for the elliptic equation (except IF sol.=0) ! (output in gcx with boundary conditions applied) kindic = 0 IF( ncut == 0 ) THEN IF( nsolv == 1 ) THEN ! diagonal preconditioned conjuguate gradient CALL sol_pcg( kindic ) ELSEIF( nsolv == 2 ) THEN ! successive-over-relaxation CALL sol_sor( kindic ) ELSEIF( nsolv == 3 ) THEN ! FETI solver CALL sol_fet( kindic ) ELSE ! e r r o r in nsolv namelist parameter IF(lwp) WRITE(numout,cform_err) IF(lwp) WRITE(numout,*) ' dyn_spg_fsc : e r r o r, nsolv = 1, 2 or 3' IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~ not = ', nsolv nstop = nstop + 1 ENDIF ENDIF ! 6. Transport divergence gradient multiplied by z2dt ! -----------------------------------------------==== DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! trend of Transport divergence gradient ztdgu = znugdt * (gcx(ji+1,jj ) - gcx(ji,jj) ) / e1u(ji,jj) ztdgv = znugdt * (gcx(ji ,jj+1) - gcx(ji,jj) ) / e2v(ji,jj) ! multiplied by z2dt #if defined key_obc ! caution : grad D = 0 along open boundaries spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj) spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj) #else spgu(ji,jj) = z2dt * ztdgu spgv(ji,jj) = z2dt * ztdgv #endif #if defined key_trddyn ! save the transport divergence gradient trends utrd2(ji,jj,2) = utrd2(ji,jj,2) + ztdgu vtrd2(ji,jj,2) = vtrd2(ji,jj,2) + ztdgv #endif END DO END DO ! 7. Add the trends multiplied by z2dt to the after velocity ! ----------------------------------------------------------- ! ( c a u t i o n : (ua,va) here are the after velocity not the ! trend, the leap-frog time stepping will not ! be done in dynnxt.F routine) DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ua(ji,jj,jk) = (ua(ji,jj,jk) + spgu(ji,jj)) * umask(ji,jj,jk) va(ji,jj,jk) = (va(ji,jj,jk) + spgv(ji,jj)) * vmask(ji,jj,jk) #if defined key_trddyn ! save the surface pressure gradient trend for diagnostics utrd(ji,jj,jk,8) = utrd(ji,jj,jk,8) + spgu(ji,jj)/z2dt vtrd(ji,jj,jk,8) = vtrd(ji,jj,jk,8) + spgv(ji,jj)/z2dt #endif END DO END DO END DO IF( l_ctl .AND. lwp ) THEN ! print sum trends (used for debugging) WRITE(numout,*) ' spg - Ua: ', SUM( ua(2:jpim1,2:jpjm1,1:jpkm1)*umask(2:jpim1,2:jpjm1,1:jpkm1) ), & & ' Va: ', SUM( va(2:jpim1,2:jpjm1,1:jpkm1)*vmask(2:jpim1,2:jpjm1,1:jpkm1) ) ENDIF ! 8. Sea surface elevation time stepping ! -------------------------------------- ! Euler (forward) time stepping, no time filter IF( neuler == 0 .AND. kt == nit000 ) THEN DO jj = 2, jpjm1 !!bug if 1, jpj the lbc_lnk call can be suppress DO ji = 1, jpi ! after free surface elevation zssha = sshb(ji,jj) + rdt * ( wn(ji,jj,1) - emp(ji,jj) * zraur ) * tmask(ji,jj,1) ! swap of arrays sshb(ji,jj) = sshn(ji,jj) sshn(ji,jj) = zssha END DO END DO ELSE ! Leap-frog time stepping and time filter DO jj = 2, jpjm1 !!bug if 1, jpj the lbc_lnk call can be suppress DO ji = 1, jpi ! after free surface elevation zssha = sshb(ji,jj) + z2dt * ( wn(ji,jj,1) - emp(ji,jj) * zraur ) * tmask(ji,jj,1) ! time filter and array swap sshb(ji,jj) = atfp * ( sshb(ji,jj) + zssha ) + atfp1 * sshn(ji,jj) sshn(ji,jj) = zssha END DO END DO ENDIF ! Boundary conditions on sshn CALL lbc_lnk( sshn, 'T', 1. ) END SUBROUTINE dyn_spg_fsc #else !!---------------------------------------------------------------------- !! Default case : Empty module No standart free surface cst volume !!---------------------------------------------------------------------- LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_fsc = .FALSE. !: free surface constant volume flag CONTAINS SUBROUTINE dyn_spg_fsc( kt, kindic ) ! Empty routine WRITE(*,*) 'dyn_spg_fsc: You should not have seen this print! error?', kt, kindic END SUBROUTINE dyn_spg_fsc #endif !!====================================================================== END MODULE dynspg_fsc