[5758] | 1 | MODULE ldfc1d_c2d |
---|
| 2 | !!====================================================================== |
---|
| 3 | !! *** MODULE ldfc1d_c2d *** |
---|
| 4 | !! Ocean physics: profile and horizontal shape of lateral eddy coefficients |
---|
| 5 | !!===================================================================== |
---|
| 6 | !! History : 3.7 ! 2013-12 (G. Madec) restructuration/simplification of aht/aeiv specification, |
---|
| 7 | !! ! add velocity dependent coefficient and optional read in file |
---|
| 8 | !!---------------------------------------------------------------------- |
---|
| 9 | |
---|
| 10 | !!---------------------------------------------------------------------- |
---|
| 11 | !! ldf_c1d : ah reduced by 1/4 on the vertical (tanh profile, inflection at 300m) |
---|
| 12 | !! ldf_c2d : ah = F(e1,e2) (laplacian or = F(e1^3,e2^3) (bilaplacian) |
---|
| 13 | !!---------------------------------------------------------------------- |
---|
| 14 | USE oce ! ocean dynamics and tracers |
---|
| 15 | USE dom_oce ! ocean space and time domain |
---|
| 16 | USE phycst ! physical constants |
---|
| 17 | ! |
---|
| 18 | USE in_out_manager ! I/O manager |
---|
| 19 | USE lib_mpp ! distribued memory computing library |
---|
| 20 | USE lbclnk ! ocean lateral boundary conditions (or mpp link) |
---|
| 21 | |
---|
| 22 | IMPLICIT NONE |
---|
| 23 | PRIVATE |
---|
| 24 | |
---|
| 25 | PUBLIC ldf_c1d ! called by ldftra and ldfdyn modules |
---|
| 26 | PUBLIC ldf_c2d ! called by ldftra and ldfdyn modules |
---|
| 27 | |
---|
| 28 | |
---|
| 29 | !! * Substitutions |
---|
| 30 | # include "vectopt_loop_substitute.h90" |
---|
| 31 | !!---------------------------------------------------------------------- |
---|
| 32 | !! NEMO/OPA 3.7 , NEMO Consortium (2015) |
---|
| 33 | !! $Id: $ |
---|
| 34 | !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) |
---|
| 35 | !!---------------------------------------------------------------------- |
---|
| 36 | CONTAINS |
---|
| 37 | |
---|
| 38 | SUBROUTINE ldf_c1d( cd_type, prat, pahs1, pahs2, pah1, pah2 ) |
---|
| 39 | !!---------------------------------------------------------------------- |
---|
| 40 | !! *** ROUTINE ldf_c1d *** |
---|
| 41 | !! |
---|
| 42 | !! ** Purpose : 1D eddy diffusivity/viscosity coefficients |
---|
| 43 | !! |
---|
| 44 | !! ** Method : 1D eddy diffusivity coefficients F( depth ) |
---|
| 45 | !! Reduction by prat from surface to bottom |
---|
| 46 | !! hyperbolic tangent profile with inflection point |
---|
| 47 | !! at zh=500m and a width of zw=200m |
---|
| 48 | !! |
---|
| 49 | !! cd_type = TRA pah1, pah2 defined at U- and V-points |
---|
| 50 | !! DYN pah1, pah2 defined at T- and F-points |
---|
| 51 | !!---------------------------------------------------------------------- |
---|
| 52 | CHARACTER(len=2) , INTENT(in ) :: cd_type ! DYNamique or TRAcers |
---|
| 53 | REAL(wp) , INTENT(in ) :: prat ! ratio surface/deep values [-] |
---|
| 54 | REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pahs1, pahs2 ! surface value of eddy coefficient [m2/s] |
---|
| 55 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pah1 , pah2 ! eddy coefficient [m2/s] |
---|
| 56 | ! |
---|
| 57 | INTEGER :: ji, jj, jk ! dummy loop indices |
---|
| 58 | REAL(wp) :: zh, zc, zdep1 ! local scalars |
---|
| 59 | REAL(wp) :: zw , zdep2 ! - - |
---|
| 60 | !!---------------------------------------------------------------------- |
---|
| 61 | |
---|
| 62 | ! initialization of the profile |
---|
| 63 | zh = 500._wp ! depth of the inflection point [m] |
---|
| 64 | zw = 1._wp / 200._wp ! width^-1 - - - [1/m] |
---|
| 65 | ! ! associated coefficient [-] |
---|
| 66 | zc = ( 1._wp - prat ) / ( 1._wp + TANH( zh * zw) ) |
---|
| 67 | ! |
---|
| 68 | ! |
---|
| 69 | SELECT CASE( cd_type ) ! point of calculation |
---|
| 70 | ! |
---|
| 71 | CASE( 'DYN' ) ! T- and F-points |
---|
| 72 | DO jk = 1, jpk ! pah1 at T-point |
---|
[6140] | 73 | pah1(:,:,jk) = pahs1(:,:) * ( prat + zc * ( 1._wp + TANH( - ( gdept_n(:,:,jk) - zh ) * zw) ) ) * tmask(:,:,jk) |
---|
[5758] | 74 | END DO |
---|
| 75 | DO jk = 1, jpk ! pah2 at F-point (zdep2 is an approximation in zps-coord.) |
---|
| 76 | DO jj = 1, jpjm1 |
---|
| 77 | DO ji = 1, fs_jpim1 |
---|
[6140] | 78 | zdep2 = ( gdept_n(ji,jj+1,jk) + gdept_n(ji+1,jj+1,jk) & |
---|
| 79 | & + gdept_n(ji,jj ,jk) + gdept_n(ji+1,jj ,jk) ) * 0.25_wp |
---|
[5758] | 80 | pah2(ji,jj,jk) = pahs2(ji,jj) * ( prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) * fmask(ji,jj,jk) |
---|
| 81 | END DO |
---|
| 82 | END DO |
---|
| 83 | END DO |
---|
| 84 | CALL lbc_lnk( pah2, 'F', 1. ) ! Lateral boundary conditions |
---|
| 85 | ! |
---|
| 86 | CASE( 'TRA' ) ! U- and V-points (zdep1 & 2 are an approximation in zps-coord.) |
---|
| 87 | DO jk = 1, jpk |
---|
| 88 | DO jj = 1, jpjm1 |
---|
| 89 | DO ji = 1, fs_jpim1 |
---|
[6140] | 90 | zdep1 = ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) * 0.5_wp |
---|
| 91 | zdep2 = ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) * 0.5_wp |
---|
[5758] | 92 | pah1(ji,jj,jk) = pahs1(ji,jj) * ( prat + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) ) ) * umask(ji,jj,jk) |
---|
| 93 | pah2(ji,jj,jk) = pahs2(ji,jj) * ( prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) * vmask(ji,jj,jk) |
---|
| 94 | END DO |
---|
| 95 | END DO |
---|
| 96 | END DO |
---|
| 97 | CALL lbc_lnk( pah1, 'U', 1. ) ! Lateral boundary conditions |
---|
| 98 | CALL lbc_lnk( pah2, 'V', 1. ) |
---|
| 99 | ! |
---|
| 100 | END SELECT |
---|
| 101 | ! |
---|
| 102 | END SUBROUTINE ldf_c1d |
---|
| 103 | |
---|
| 104 | |
---|
| 105 | SUBROUTINE ldf_c2d( cd_type, cd_op, pah0, pah1, pah2 ) |
---|
| 106 | !!---------------------------------------------------------------------- |
---|
| 107 | !! *** ROUTINE ldf_c2d *** |
---|
| 108 | !! |
---|
| 109 | !! ** Purpose : 2D eddy diffusivity/viscosity coefficients |
---|
| 110 | !! |
---|
| 111 | !! ** Method : 2D eddy diffusivity coefficients F( e1 , e2 ) |
---|
| 112 | !! laplacian operator : ah proportional to the scale factor [m2/s] |
---|
| 113 | !! bilaplacian operator : ah proportional to the (scale factor)^3 [m4/s] |
---|
| 114 | !! In both cases, pah0 is the maximum value reached by the coefficient |
---|
| 115 | !! at the Equator in case of e1=ra*rad= ~111km, not over the whole domain. |
---|
| 116 | !! |
---|
| 117 | !! cd_type = TRA pah1, pah2 defined at U- and V-points |
---|
| 118 | !! DYN pah1, pah2 defined at T- and F-points |
---|
| 119 | !!---------------------------------------------------------------------- |
---|
| 120 | CHARACTER(len=3) , INTENT(in ) :: cd_type ! DYNamique or TRAcers |
---|
| 121 | CHARACTER(len=3) , INTENT(in ) :: cd_op ! operator: LAPlacian BiLaPlacian |
---|
| 122 | REAL(wp) , INTENT(in ) :: pah0 ! eddy coefficient [m2/s or m4/s] |
---|
| 123 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: pah1, pah2 ! eddy coefficient [m2/s or m4/s] |
---|
| 124 | ! |
---|
| 125 | INTEGER :: ji, jj, jk ! dummy loop indices |
---|
| 126 | REAL(wp) :: za00, zd_max, zemax1, zemax2 ! local scalar |
---|
| 127 | !!---------------------------------------------------------------------- |
---|
| 128 | ! |
---|
| 129 | zd_max = ra * rad ! = 1 degree at the equator in meters |
---|
| 130 | ! |
---|
| 131 | IF(lwp) THEN |
---|
| 132 | WRITE(numout,*) 'ldf_c2d : aht = rn_aht0 * max(e1,e2)/e_equator ( laplacian) ' |
---|
| 133 | WRITE(numout,*) '~~~~~~~ or = rn_bht0 * [max(e1,e2)/e_equator]**3 (bilaplacian)' |
---|
| 134 | WRITE(numout,*) |
---|
| 135 | ENDIF |
---|
| 136 | ! |
---|
| 137 | SELECT CASE( cd_type ) !== surface values ==! (depending on DYN/TRA) |
---|
| 138 | ! |
---|
| 139 | CASE( 'DYN' ) ! T- and F-points |
---|
| 140 | IF( cd_op == 'LAP' ) THEN ! laplacian operator |
---|
| 141 | IF(lwp) WRITE(numout,*) ' momentum laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)' |
---|
| 142 | za00 = pah0 / zd_max |
---|
| 143 | DO jj = 1, jpj |
---|
| 144 | DO ji = 1, jpi |
---|
| 145 | zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1) |
---|
| 146 | zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1) |
---|
| 147 | pah1(ji,jj,1) = za00 * zemax1 |
---|
| 148 | pah2(ji,jj,1) = za00 * zemax2 |
---|
| 149 | END DO |
---|
| 150 | END DO |
---|
| 151 | ELSEIF( cd_op == 'BLP' ) THEN ! bilaplacian operator |
---|
| 152 | IF(lwp) WRITE(numout,*) ' momentum bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3' |
---|
| 153 | za00 = pah0 / ( zd_max * zd_max * zd_max ) |
---|
| 154 | DO jj = 1, jpj |
---|
| 155 | DO ji = 1, jpi |
---|
| 156 | zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1) |
---|
| 157 | zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1) |
---|
| 158 | pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1 |
---|
| 159 | pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2 |
---|
| 160 | END DO |
---|
| 161 | END DO |
---|
| 162 | ELSE ! NO diffusion/viscosity |
---|
| 163 | CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' ) |
---|
| 164 | ENDIF |
---|
| 165 | ! ! deeper values (LAP and BLP cases) |
---|
| 166 | DO jk = 2, jpk |
---|
| 167 | pah1(:,:,jk) = pah1(:,:,1) * tmask(:,:,jk) |
---|
| 168 | pah2(:,:,jk) = pah2(:,:,1) * fmask(:,:,jk) |
---|
| 169 | END DO |
---|
| 170 | ! |
---|
| 171 | CASE( 'TRA' ) ! U- and V-points (approximation in zps-coord.) |
---|
| 172 | IF( cd_op == 'LAP' ) THEN ! laplacian operator |
---|
| 173 | IF(lwp) WRITE(numout,*) ' tracer laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)' |
---|
| 174 | za00 = pah0 / zd_max |
---|
| 175 | DO jj = 1, jpj |
---|
| 176 | DO ji = 1, jpi |
---|
| 177 | zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1) |
---|
| 178 | zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1) |
---|
| 179 | pah1(ji,jj,1) = za00 * zemax1 |
---|
| 180 | pah2(ji,jj,1) = za00 * zemax2 |
---|
| 181 | END DO |
---|
| 182 | END DO |
---|
| 183 | ELSEIF( cd_op == 'BLP' ) THEN ! bilaplacian operator (NB: square root of the coeff) |
---|
| 184 | IF(lwp) WRITE(numout,*) ' tracer bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3' |
---|
| 185 | za00 = pah0 / ( zd_max * zd_max * zd_max ) |
---|
| 186 | DO jj = 1, jpj |
---|
| 187 | DO ji = 1, jpi |
---|
| 188 | zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1) |
---|
| 189 | zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1) |
---|
| 190 | pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1 |
---|
| 191 | pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2 |
---|
| 192 | END DO |
---|
| 193 | END DO |
---|
| 194 | ELSE ! NO diffusion/viscosity |
---|
| 195 | CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' ) |
---|
| 196 | ENDIF |
---|
| 197 | ! ! deeper values (LAP and BLP cases) |
---|
| 198 | DO jk = 2, jpk |
---|
| 199 | pah1(:,:,jk) = pah1(:,:,1) * umask(:,:,jk) |
---|
| 200 | pah2(:,:,jk) = pah2(:,:,1) * vmask(:,:,jk) |
---|
| 201 | END DO |
---|
| 202 | ! |
---|
| 203 | END SELECT |
---|
| 204 | ! |
---|
| 205 | END SUBROUTINE ldf_c2d |
---|
| 206 | |
---|
| 207 | !!====================================================================== |
---|
| 208 | END MODULE ldfc1d_c2d |
---|