[5777] | 1 | MODULE dynldf_lap_blp |
---|
[3] | 2 | !!====================================================================== |
---|
[5777] | 3 | !! *** MODULE dynldf_lap_blp *** |
---|
| 4 | !! Ocean dynamics: lateral viscosity trend (laplacian and bilaplacian) |
---|
[3] | 5 | !!====================================================================== |
---|
[6140] | 6 | !! History : 3.7 ! 2014-01 (G. Madec, S. Masson) Original code, re-entrant laplacian |
---|
[2715] | 7 | !!---------------------------------------------------------------------- |
---|
[3] | 8 | |
---|
| 9 | !!---------------------------------------------------------------------- |
---|
[5777] | 10 | !! dyn_ldf_lap : update the momentum trend with the lateral viscosity using an iso-level laplacian operator |
---|
| 11 | !! dyn_ldf_blp : update the momentum trend with the lateral viscosity using an iso-level bilaplacian operator |
---|
[3] | 12 | !!---------------------------------------------------------------------- |
---|
[5777] | 13 | USE oce ! ocean dynamics and tracers |
---|
| 14 | USE dom_oce ! ocean space and time domain |
---|
| 15 | USE ldfdyn ! lateral diffusion: eddy viscosity coef. |
---|
| 16 | USE ldfslp ! iso-neutral slopes |
---|
| 17 | USE zdf_oce ! ocean vertical physics |
---|
[4990] | 18 | ! |
---|
[5777] | 19 | USE in_out_manager ! I/O manager |
---|
| 20 | USE lbclnk ! ocean lateral boundary conditions (or mpp link) |
---|
| 21 | USE timing ! Timing |
---|
[3] | 22 | |
---|
| 23 | IMPLICIT NONE |
---|
| 24 | PRIVATE |
---|
| 25 | |
---|
[5777] | 26 | PUBLIC dyn_ldf_lap ! called by dynldf.F90 |
---|
| 27 | PUBLIC dyn_ldf_blp ! called by dynldf.F90 |
---|
[3] | 28 | |
---|
| 29 | !! * Substitutions |
---|
| 30 | # include "vectopt_loop_substitute.h90" |
---|
| 31 | !!---------------------------------------------------------------------- |
---|
[5777] | 32 | !! NEMO/OPA 3.7 , NEMO Consortium (2014) |
---|
[1152] | 33 | !! $Id$ |
---|
[2715] | 34 | !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) |
---|
[3] | 35 | !!---------------------------------------------------------------------- |
---|
| 36 | CONTAINS |
---|
| 37 | |
---|
[5777] | 38 | SUBROUTINE dyn_ldf_lap( kt, pub, pvb, pua, pva, kpass ) |
---|
[3] | 39 | !!---------------------------------------------------------------------- |
---|
| 40 | !! *** ROUTINE dyn_ldf_lap *** |
---|
| 41 | !! |
---|
[5777] | 42 | !! ** Purpose : Compute the before horizontal momentum diffusive |
---|
| 43 | !! trend and add it to the general trend of momentum equation. |
---|
[3] | 44 | !! |
---|
[5777] | 45 | !! ** Method : The Laplacian operator apply on horizontal velocity is |
---|
| 46 | !! writen as : grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) |
---|
[3] | 47 | !! |
---|
[5777] | 48 | !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb. |
---|
[3] | 49 | !!---------------------------------------------------------------------- |
---|
[5777] | 50 | INTEGER , INTENT(in ) :: kt ! ocean time-step index |
---|
| 51 | INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage |
---|
| 52 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pub, pvb ! before velocity [m/s] |
---|
| 53 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! velocity trend [m/s2] |
---|
[2715] | 54 | ! |
---|
[5777] | 55 | INTEGER :: ji, jj, jk ! dummy loop indices |
---|
| 56 | REAL(wp) :: zsign ! local scalars |
---|
| 57 | REAL(wp) :: zua, zva ! local scalars |
---|
[7910] | 58 | REAL(wp), DIMENSION(jpi,jpj) :: zcur, zdiv |
---|
[3] | 59 | !!---------------------------------------------------------------------- |
---|
[2715] | 60 | ! |
---|
[5777] | 61 | IF( kt == nit000 .AND. lwp ) THEN |
---|
| 62 | WRITE(numout,*) |
---|
| 63 | WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass |
---|
| 64 | WRITE(numout,*) '~~~~~~~ ' |
---|
| 65 | ENDIF |
---|
[3294] | 66 | ! |
---|
[5777] | 67 | IF( nn_timing == 1 ) CALL timing_start('dyn_ldf_lap') |
---|
| 68 | ! |
---|
| 69 | ! |
---|
| 70 | IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign |
---|
| 71 | ELSE ; zsign = -1._wp ! (eddy viscosity coef. >0) |
---|
[3] | 72 | ENDIF |
---|
[5777] | 73 | ! |
---|
[3] | 74 | ! ! =============== |
---|
| 75 | DO jk = 1, jpkm1 ! Horizontal slab |
---|
| 76 | ! ! =============== |
---|
[5777] | 77 | DO jj = 2, jpj |
---|
| 78 | DO ji = fs_2, jpi ! vector opt. |
---|
| 79 | ! ! ahm * e3 * curl (computed from 1 to jpim1/jpjm1) |
---|
[6140] | 80 | !!gm open question here : e3f at before or now ? probably now... |
---|
[5860] | 81 | !!gm note that ahmf has already been multiplied by fmask |
---|
[6140] | 82 | zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & |
---|
[5860] | 83 | & * ( e2v(ji ,jj-1) * pvb(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk) & |
---|
[6140] | 84 | & - e1u(ji-1,jj ) * pub(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk) ) |
---|
[5777] | 85 | ! ! ahm * div (computed from 2 to jpi/jpj) |
---|
[6140] | 86 | !!gm note that ahmt has already been multiplied by tmask |
---|
| 87 | zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk) & |
---|
| 88 | & * ( e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk) & |
---|
| 89 | & + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk) ) |
---|
[5777] | 90 | END DO |
---|
| 91 | END DO |
---|
| 92 | ! |
---|
| 93 | DO jj = 2, jpjm1 ! - curl( curl) + grad( div ) |
---|
[3] | 94 | DO ji = fs_2, fs_jpim1 ! vector opt. |
---|
[6140] | 95 | pua(ji,jj,jk) = pua(ji,jj,jk) + zsign * ( & |
---|
| 96 | & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk) & |
---|
[5777] | 97 | & + ( zdiv(ji+1,jj) - zdiv(ji,jj ) ) * r1_e1u(ji,jj) ) |
---|
| 98 | ! |
---|
[6140] | 99 | pva(ji,jj,jk) = pva(ji,jj,jk) + zsign * ( & |
---|
| 100 | & ( zcur(ji,jj ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk) & |
---|
[5777] | 101 | & + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) ) |
---|
[3] | 102 | END DO |
---|
| 103 | END DO |
---|
| 104 | ! ! =============== |
---|
| 105 | END DO ! End of slab |
---|
| 106 | ! ! =============== |
---|
[5777] | 107 | ! |
---|
[3294] | 108 | IF( nn_timing == 1 ) CALL timing_stop('dyn_ldf_lap') |
---|
| 109 | ! |
---|
[3] | 110 | END SUBROUTINE dyn_ldf_lap |
---|
| 111 | |
---|
[5777] | 112 | |
---|
| 113 | SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva ) |
---|
| 114 | !!---------------------------------------------------------------------- |
---|
| 115 | !! *** ROUTINE dyn_ldf_blp *** |
---|
| 116 | !! |
---|
| 117 | !! ** Purpose : Compute the before lateral momentum viscous trend |
---|
| 118 | !! and add it to the general trend of momentum equation. |
---|
| 119 | !! |
---|
| 120 | !! ** Method : The lateral viscous trends is provided by a bilaplacian |
---|
| 121 | !! operator applied to before field (forward in time). |
---|
| 122 | !! It is computed by two successive calls to dyn_ldf_lap routine |
---|
| 123 | !! |
---|
| 124 | !! ** Action : pta updated with the before rotated bilaplacian diffusion |
---|
| 125 | !!---------------------------------------------------------------------- |
---|
| 126 | INTEGER , INTENT(in ) :: kt ! ocean time-step index |
---|
| 127 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pub, pvb ! before velocity fields |
---|
| 128 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! momentum trend |
---|
| 129 | ! |
---|
[7910] | 130 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zulap, zvlap ! laplacian at u- and v-point |
---|
[5777] | 131 | !!---------------------------------------------------------------------- |
---|
| 132 | ! |
---|
| 133 | IF( nn_timing == 1 ) CALL timing_start('dyn_ldf_blp') |
---|
| 134 | ! |
---|
| 135 | ! |
---|
| 136 | IF( kt == nit000 ) THEN |
---|
| 137 | IF(lwp) WRITE(numout,*) |
---|
| 138 | IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum ' |
---|
| 139 | IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' |
---|
| 140 | ENDIF |
---|
| 141 | ! |
---|
[7753] | 142 | zulap(:,:,:) = 0._wp |
---|
| 143 | zvlap(:,:,:) = 0._wp |
---|
[5777] | 144 | ! |
---|
| 145 | CALL dyn_ldf_lap( kt, pub, pvb, zulap, zvlap, 1 ) ! rotated laplacian applied to ptb (output in zlap) |
---|
| 146 | ! |
---|
| 147 | CALL lbc_lnk( zulap(:,:,:) , 'U', -1. ) ! Lateral boundary conditions |
---|
| 148 | CALL lbc_lnk( zvlap(:,:,:) , 'V', -1. ) |
---|
| 149 | ! |
---|
| 150 | CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 ) ! rotated laplacian applied to zlap (output in pta) |
---|
| 151 | ! |
---|
| 152 | ! |
---|
| 153 | IF( nn_timing == 1 ) CALL timing_stop('dyn_ldf_blp') |
---|
| 154 | ! |
---|
| 155 | END SUBROUTINE dyn_ldf_blp |
---|
| 156 | |
---|
[3] | 157 | !!====================================================================== |
---|
[5777] | 158 | END MODULE dynldf_lap_blp |
---|