Changeset 14072 for NEMO/trunk/src/OCE/TRA
- Timestamp:
- 2020-12-04T08:48:38+01:00 (3 years ago)
- Location:
- NEMO/trunk/src/OCE/TRA
- Files:
-
- 22 edited
- 2 copied
-
eosbn2.F90 (modified) (19 diffs)
-
traadv.F90 (modified) (6 diffs)
-
traadv_cen.F90 (modified) (7 diffs)
-
traadv_fct.F90 (modified) (25 diffs)
-
traadv_fct_lf.F90 (copied) (copied from NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/TRA/traadv_fct_lf.F90)
-
traadv_mus.F90 (modified) (6 diffs)
-
traadv_mus_lf.F90 (copied) (copied from NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/TRA/traadv_mus_lf.F90)
-
traadv_qck.F90 (modified) (15 diffs)
-
traadv_ubs.F90 (modified) (14 diffs)
-
traatf.F90 (modified) (19 diffs)
-
traatf_qco.F90 (modified) (1 diff)
-
trabbc.F90 (modified) (5 diffs)
-
trabbl.F90 (modified) (4 diffs)
-
tradmp.F90 (modified) (5 diffs)
-
traisf.F90 (modified) (6 diffs)
-
traldf.F90 (modified) (4 diffs)
-
traldf_iso.F90 (modified) (12 diffs)
-
traldf_lap_blp.F90 (modified) (8 diffs)
-
traldf_triad.F90 (modified) (5 diffs)
-
tranpc.F90 (modified) (10 diffs)
-
traqsr.F90 (modified) (14 diffs)
-
trasbc.F90 (modified) (9 diffs)
-
trazdf.F90 (modified) (7 diffs)
-
zpshde.F90 (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/TRA/eosbn2.F90
r14010 r14072 31 31 !! bn2 : compute the Brunt-Vaisala frequency 32 32 !! eos_pt_from_ct: compute the potential temperature from the Conservative Temperature 33 !! eos_rab : generic interface of in situ thermal/haline expansion ratio 33 !! eos_rab : generic interface of in situ thermal/haline expansion ratio 34 34 !! eos_rab_3d : compute in situ thermal/haline expansion ratio 35 35 !! eos_rab_2d : compute in situ thermal/haline expansion ratio for 2d fields … … 46 46 USE in_out_manager ! I/O manager 47 47 USE lib_mpp ! MPP library 48 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 48 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 49 49 USE prtctl ! Print control 50 50 USE lbclnk ! ocean lateral boundary conditions … … 63 63 END INTERFACE 64 64 ! 65 INTERFACE eos_fzp 65 INTERFACE eos_fzp 66 66 MODULE PROCEDURE eos_fzp_2d, eos_fzp_0d 67 67 END INTERFACE … … 89 89 90 90 ! !!! simplified eos coefficients (default value: Vallis 2006) 91 REAL(wp) :: rn_a0 = 1.6550e-1_wp ! thermal expansion coeff. 92 REAL(wp) :: rn_b0 = 7.6554e-1_wp ! saline expansion coeff. 93 REAL(wp) :: rn_lambda1 = 5.9520e-2_wp ! cabbeling coeff. in T^2 94 REAL(wp) :: rn_lambda2 = 5.4914e-4_wp ! cabbeling coeff. in S^2 95 REAL(wp) :: rn_mu1 = 1.4970e-4_wp ! thermobaric coeff. in T 96 REAL(wp) :: rn_mu2 = 1.1090e-5_wp ! thermobaric coeff. in S 97 REAL(wp) :: rn_nu = 2.4341e-3_wp ! cabbeling coeff. in theta*salt 98 91 REAL(wp) :: rn_a0 = 1.6550e-1_wp ! thermal expansion coeff. 92 REAL(wp) :: rn_b0 = 7.6554e-1_wp ! saline expansion coeff. 93 REAL(wp) :: rn_lambda1 = 5.9520e-2_wp ! cabbeling coeff. in T^2 94 REAL(wp) :: rn_lambda2 = 5.4914e-4_wp ! cabbeling coeff. in S^2 95 REAL(wp) :: rn_mu1 = 1.4970e-4_wp ! thermobaric coeff. in T 96 REAL(wp) :: rn_mu2 = 1.1090e-5_wp ! thermobaric coeff. in S 97 REAL(wp) :: rn_nu = 2.4341e-3_wp ! cabbeling coeff. in theta*salt 98 99 99 ! TEOS10/EOS80 parameters 100 100 REAL(wp) :: r1_S0, r1_T0, r1_Z0, rdeltaS 101 101 102 102 ! EOS parameters 103 103 REAL(wp) :: EOS000 , EOS100 , EOS200 , EOS300 , EOS400 , EOS500 , EOS600 … … 117 117 REAL(wp) :: EOS022 118 118 REAL(wp) :: EOS003 , EOS103 119 REAL(wp) :: EOS013 120 119 REAL(wp) :: EOS013 120 121 121 ! ALPHA parameters 122 122 REAL(wp) :: ALP000 , ALP100 , ALP200 , ALP300 , ALP400 , ALP500 … … 133 133 REAL(wp) :: ALP012 134 134 REAL(wp) :: ALP003 135 135 136 136 ! BETA parameters 137 137 REAL(wp) :: BET000 , BET100 , BET200 , BET300 , BET400 , BET500 … … 160 160 REAL(wp) :: PEN002 , PEN102 161 161 REAL(wp) :: PEN012 162 162 163 163 ! ALPHA_PEN parameters 164 164 REAL(wp) :: APE000 , APE100 , APE200 , APE300 … … 295 295 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs & 296 296 & - rn_nu * zt * zs 297 ! 297 ! 298 298 prd(ji,jj,jk) = zn * r1_rho0 * ztm ! density anomaly (masked) 299 299 END_3D … … 448 448 END_3D 449 449 ENDIF 450 450 451 451 CASE( np_seos ) !== simplified EOS ==! 452 452 ! … … 997 997 !! *** ROUTINE bn2 *** 998 998 !! 999 !! ** Purpose : Compute the local Brunt-Vaisala frequency at the 999 !! ** Purpose : Compute the local Brunt-Vaisala frequency at the 1000 1000 !! time-step of the input arguments 1001 1001 !! … … 1004 1004 !! N.B. N^2 is set one for all to zero at jk=1 in istate module. 1005 1005 !! 1006 !! ** Action : pn2 : square of the brunt-vaisala frequency at w-point 1006 !! ** Action : pn2 : square of the brunt-vaisala frequency at w-point 1007 1007 !! 1008 1008 !!---------------------------------------------------------------------- … … 1021 1021 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1 ) ! interior points only (2=< jk =< jpkm1 ); surface and bottom value set to zero one for all in istate.F90 1022 1022 zrw = ( gdepw(ji,jj,jk ,Kmm) - gdept(ji,jj,jk,Kmm) ) & 1023 & / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 1024 ! 1025 zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 1023 & / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 1024 ! 1025 zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 1026 1026 zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw 1027 1027 ! … … 1151 1151 CALL ctl_stop( 'eos_fzp_2d:', ctmp1 ) 1152 1152 ! 1153 END SELECT 1153 END SELECT 1154 1154 ! 1155 1155 END SUBROUTINE eos_fzp_2d_t … … 1208 1208 !! ** Purpose : Calculates nonlinear anomalies of alpha_PE, beta_PE and PE at T-points 1209 1209 !! 1210 !! ** Method : PE is defined analytically as the vertical 1210 !! ** Method : PE is defined analytically as the vertical 1211 1211 !! primitive of EOS times -g integrated between 0 and z>0. 1212 1212 !! pen is the nonlinear bsq-PE anomaly: pen = ( PE - rho0 gz ) / rho0 gz - rd 1213 !! = 1/z * /int_0^z rd dz - rd 1213 !! = 1/z * /int_0^z rd dz - rd 1214 1214 !! where rd is the density anomaly (see eos_rhd function) 1215 1215 !! ab_pe are partial derivatives of PE anomaly with respect to T and S: … … 1275 1275 ! 1276 1276 zn = ( zn2 * zh + zn1 ) * zh + zn0 1277 ! 1277 ! 1278 1278 pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm 1279 1279 ! … … 1290 1290 ! 1291 1291 zn = ( zn2 * zh + zn1 ) * zh + zn0 1292 ! 1292 ! 1293 1293 pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm 1294 1294 ! … … 1370 1370 IF(lwp) WRITE(numout,*) ' ==>>> use of TEOS-10 equation of state (cons. temp. and abs. salinity)' 1371 1371 ! 1372 l_useCT = .TRUE. ! model temperature is Conservative temperature 1372 l_useCT = .TRUE. ! model temperature is Conservative temperature 1373 1373 ! 1374 1374 rdeltaS = 32._wp … … 1751 1751 1752 1752 r1_S0 = 0.875_wp/35.16504_wp ! Used to convert CT in potential temperature when using bulk formulae (eos_pt_from_ct) 1753 1753 1754 1754 IF(lwp) THEN 1755 1755 WRITE(numout,*) … … 1775 1775 END SELECT 1776 1776 ! 1777 rho0_rcp = rho0 * rcp 1777 rho0_rcp = rho0 * rcp 1778 1778 r1_rho0 = 1._wp / rho0 1779 1779 r1_rcp = 1._wp / rcp 1780 r1_rho0_rcp = 1._wp / rho0_rcp 1780 r1_rho0_rcp = 1._wp / rho0_rcp 1781 1781 ! 1782 1782 IF(lwp) THEN -
NEMO/trunk/src/OCE/TRA/traadv.F90
r13982 r14072 2 2 !!============================================================================== 3 3 !! *** MODULE traadv *** 4 !! Ocean active tracers: advection trend 4 !! Ocean active tracers: advection trend 5 5 !!============================================================================== 6 6 !! History : 2.0 ! 2005-11 (G. Madec) Original code 7 7 !! 3.3 ! 2010-09 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport 8 8 !! 3.6 ! 2011-06 (G. Madec) Addition of Mixed Layer Eddy parameterisation 9 !! 3.7 ! 2014-05 (G. Madec) Add 2nd/4th order cases for CEN and FCT schemes 9 !! 3.7 ! 2014-05 (G. Madec) Add 2nd/4th order cases for CEN and FCT schemes 10 10 !! - ! 2014-12 (G. Madec) suppression of cross land advection option 11 11 !! 3.6 ! 2015-06 (E. Clementi) Addition of Stokes drift in case of wave coupling … … 34 34 USE ldfslp ! Lateral diffusion: slopes of neutral surfaces 35 35 USE trd_oce ! trends: ocean variables 36 USE trdtra ! trends manager: tracers 37 USE diaptr ! Poleward heat transport 36 USE trdtra ! trends manager: tracers 37 USE diaptr ! Poleward heat transport 38 38 ! 39 39 USE in_out_manager ! I/O manager … … 195 195 CASE ( np_MUS ) ! MUSCL 196 196 ! NOTE: [tiling-comms-merge] I added this lbc_lnk as it did not validate against the trunk when using ln_zco 197 IF (nn_hls.EQ.2) THEN 197 IF (nn_hls.EQ.2) THEN 198 198 CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.) 199 199 #if defined key_loop_fusion 200 CALL tra_adv_mus_lf ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 200 CALL tra_adv_mus_lf ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 201 201 #else 202 CALL tra_adv_mus ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 202 CALL tra_adv_mus ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 203 203 #endif 204 204 ELSE 205 CALL tra_adv_mus ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 205 CALL tra_adv_mus ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 206 206 END IF 207 207 CASE ( np_UBS ) ! UBS … … 248 248 !!--------------------------------------------------------------------- 249 249 !! *** ROUTINE tra_adv_init *** 250 !! 251 !! ** Purpose : Control the consistency between namelist options for 250 !! 251 !! ** Purpose : Control the consistency between namelist options for 252 252 !! tracer advection schemes and set nadv 253 253 !!---------------------------------------------------------------------- … … 290 290 ! 291 291 ! !== Parameter control & set nadv ==! 292 ioptio = 0 292 ioptio = 0 293 293 IF( ln_traadv_OFF ) THEN ; ioptio = ioptio + 1 ; nadv = np_NO_adv ; ENDIF 294 294 IF( ln_traadv_cen ) THEN ; ioptio = ioptio + 1 ; nadv = np_CEN ; ENDIF … … 319 319 ENDIF 320 320 ! 321 ! !== Print the choice ==! 321 ! !== Print the choice ==! 322 322 IF(lwp) THEN 323 323 WRITE(numout,*) -
NEMO/trunk/src/OCE/TRA/traadv_cen.F90
r13982 r14072 13 13 USE dom_oce ! ocean space and time domain 14 14 USE eosbn2 ! equation of state 15 USE traadv_fct ! acces to routine interp_4th_cpt 15 USE traadv_fct ! acces to routine interp_4th_cpt 16 16 USE trd_oce ! trends: ocean variables 17 USE trdtra ! trends manager: tracers 17 USE trdtra ! trends manager: tracers 18 18 USE diaptr ! poleward transport diagnostics 19 19 USE diaar5 ! AR5 diagnostics … … 28 28 29 29 PUBLIC tra_adv_cen ! called by traadv.F90 30 30 31 31 REAL(wp) :: r1_6 = 1._wp / 6._wp ! =1/6 32 32 … … 46 46 47 47 SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pU, pV, pW, & 48 & Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v ) 48 & Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v ) 49 49 !!---------------------------------------------------------------------- 50 50 !! *** ROUTINE tra_adv_cen *** 51 !! 51 !! 52 52 !! ** Purpose : Compute the now trend due to the advection of tracers 53 53 !! and add it to the general trend of passive tracer equations. 54 54 !! 55 55 !! ** Method : The advection is evaluated by a 2nd or 4th order scheme 56 !! using now fields (leap-frog scheme). 56 !! using now fields (leap-frog scheme). 57 57 !! kn_cen_h = 2 ==>> 2nd order centered scheme on the horizontal 58 58 !! = 4 ==>> 4th order - - - - … … 98 98 ENDIF 99 99 ! 100 ! 100 ! 101 101 zwz(:,:, 1 ) = 0._wp ! surface & bottom vertical flux set to zero for all tracers 102 102 zwz(:,:,jpk) = 0._wp … … 155 155 IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean) 156 156 DO_2D( 1, 1, 1, 1 ) 157 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) 157 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) 158 158 END_2D 159 159 ELSE ! no ice-shelf cavities (only ocean surface) … … 163 163 ENDIF 164 164 ENDIF 165 ! 165 ! 166 166 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Divergence of advective fluxes --! 167 167 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) & … … 185 185 ! 186 186 END SUBROUTINE tra_adv_cen 187 187 188 188 !!====================================================================== 189 189 END MODULE traadv_cen -
NEMO/trunk/src/OCE/TRA/traadv_fct.F90
r13982 r14072 10 10 !! tra_adv_fct : update the tracer trend with a 3D advective trends using a 2nd or 4th order FCT scheme 11 11 !! with sub-time-stepping in the vertical direction 12 !! nonosc : compute monotonic tracer fluxes by a non-oscillatory algorithm 12 !! nonosc : compute monotonic tracer fluxes by a non-oscillatory algorithm 13 13 !! interp_4th_cpt : 4th order compact scheme for the vertical component of the advection 14 14 !!---------------------------------------------------------------------- … … 24 24 ! 25 25 USE in_out_manager ! I/O manager 26 USE iom ! 26 USE iom ! 27 27 USE lib_mpp ! MPP library 28 USE lbclnk ! ocean lateral boundary condition (or mpp link) 29 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 28 USE lbclnk ! ocean lateral boundary condition (or mpp link) 29 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 30 30 31 31 IMPLICIT NONE … … 60 60 !!---------------------------------------------------------------------- 61 61 !! *** ROUTINE tra_adv_fct *** 62 !! 62 !! 63 63 !! ** Purpose : Compute the now trend due to total advection of tracers 64 64 !! and add it to the general trend of tracer equations … … 66 66 !! ** Method : - 2nd or 4th FCT scheme on the horizontal direction 67 67 !! (choice through the value of kn_fct) 68 !! - on the vertical the 4th order is a compact scheme 69 !! - corrected flux (monotonic correction) 68 !! - on the vertical the 4th order is a compact scheme 69 !! - corrected flux (monotonic correction) 70 70 !! 71 71 !! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends … … 154 154 ! 155 155 ! !== upstream advection with initial mass fluxes & intermediate update ==! 156 ! !* upstream tracer flux in the i and j direction 156 ! !* upstream tracer flux in the i and j direction 157 157 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 158 158 ! upstream scheme … … 173 173 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface 174 174 DO_2D( 1, 1, 1, 1 ) 175 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface 175 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface 176 176 END_2D 177 177 ELSE ! no cavities: only at the ocean surface … … 181 181 ENDIF 182 182 ENDIF 183 ! 183 ! 184 184 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* trend and after field with monotonic scheme 185 185 ! ! total intermediate advective trends … … 193 193 & / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 194 194 END_3D 195 195 196 196 IF ( ll_zAimp ) THEN 197 197 CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 ) … … 215 215 END IF 216 216 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 217 IF( l_ptr ) zptry(:,:,:) = zwy(:,:,:) 217 IF( l_ptr ) zptry(:,:,:) = zwy(:,:,:) 218 218 ! 219 219 ! !== anti-diffusive flux : high order minus low order ==! … … 268 268 zC4t_u = zC2t_u + r1_6 * ( ztu(ji-1,jj ,jk) - ztu(ji+1,jj ,jk) ) 269 269 zC4t_v = zC2t_v + r1_6 * ( ztv(ji ,jj-1,jk) - ztv(ji ,jj+1,jk) ) 270 ! ! C4 minus upstream advective fluxes 270 ! ! C4 minus upstream advective fluxes 271 271 zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk) 272 272 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) … … 275 275 ! 276 276 END SELECT 277 ! 277 ! 278 278 SELECT CASE( kn_fct_v ) !* vertical anti-diffusive fluxes (w-masked interior values) 279 279 ! … … 384 384 DEALLOCATE( ztrdx, ztrdy, ztrdz ) 385 385 ENDIF 386 IF( l_ptr ) THEN 386 IF( l_ptr ) THEN 387 387 DEALLOCATE( zptry ) 388 388 ENDIF … … 394 394 !!--------------------------------------------------------------------- 395 395 !! *** ROUTINE nonosc *** 396 !! 397 !! ** Purpose : compute monotonic tracer fluxes from the upstream 398 !! scheme and the before field by a nonoscillatory algorithm 396 !! 397 !! ** Purpose : compute monotonic tracer fluxes from the upstream 398 !! scheme and the before field by a nonoscillatory algorithm 399 399 !! 400 400 !! ** Method : ... ??? … … 492 492 !!---------------------------------------------------------------------- 493 493 !! *** ROUTINE interp_4th_cpt_org *** 494 !! 494 !! 495 495 !! ** Purpose : Compute the interpolation of tracer at w-point 496 496 !! … … 503 503 REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt 504 504 !!---------------------------------------------------------------------- 505 505 506 506 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) !== build the three diagonal matrix ==! 507 507 zwd (ji,jj,jk) = 4._wp … … 514 514 zwi (ji,jj,jk) = 0._wp 515 515 zws (ji,jj,jk) = 0._wp 516 zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 516 zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 517 517 ENDIF 518 518 END_3D … … 538 538 END_2D 539 539 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 540 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 540 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 541 541 END_3D 542 542 … … 547 547 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 548 548 END_3D 549 ! 549 ! 550 550 END SUBROUTINE interp_4th_cpt_org 551 551 552 552 553 553 SUBROUTINE interp_4th_cpt( pt_in, pt_out ) 554 554 !!---------------------------------------------------------------------- 555 555 !! *** ROUTINE interp_4th_cpt *** 556 !! 556 !! 557 557 !! ** Purpose : Compute the interpolation of tracer at w-point 558 558 !! … … 582 582 ! CASE( np_CEN2 ) ! 2nd order centered at top & bottom 583 583 ! END SELECT 584 !!gm 584 !!gm 585 585 ! 586 586 IF ( ln_isfcav ) THEN ! set level two values which may not be set in ISF case … … 600 600 zwi (ji,jj,ikb) = 0._wp 601 601 zws (ji,jj,ikb) = 0._wp 602 zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) ) 602 zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) ) 603 603 END_2D 604 604 ! … … 616 616 END_2D 617 617 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) 618 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 618 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 619 619 END_3D 620 620 … … 625 625 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 626 626 END_3D 627 ! 627 ! 628 628 END SUBROUTINE interp_4th_cpt 629 629 … … 632 632 !!---------------------------------------------------------------------- 633 633 !! *** ROUTINE tridia_solver *** 634 !! 634 !! 635 635 !! ** Purpose : solve a symmetric 3diagonal system 636 636 !! 637 637 !! ** Method : solve M.t_out = RHS(t) where M is a tri diagonal matrix ( jpk*jpk ) 638 !! 638 !! 639 639 !! ( D_1 U_1 0 0 0 )( t_1 ) ( RHS_1 ) 640 640 !! ( L_2 D_2 U_2 0 0 )( t_2 ) ( RHS_2 ) … … 642 642 !! ( ... )( ... ) ( ... ) 643 643 !! ( 0 0 0 L_k D_k )( t_k ) ( RHS_k ) 644 !! 644 !! 645 645 !! M is decomposed in the product of an upper and lower triangular matrix. 646 !! The tri-diagonals matrix is given as input 3D arrays: pD, pU, pL 646 !! The tri-diagonals matrix is given as input 3D arrays: pD, pU, pL 647 647 !! (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal). 648 648 !! The solution is pta. … … 672 672 END_2D 673 673 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, kstart+1, jpkm1 ) 674 pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 674 pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 675 675 END_3D 676 676 -
NEMO/trunk/src/OCE/TRA/traadv_mus.F90
r13982 r14072 29 29 USE in_out_manager ! I/O manager 30 30 USE lib_mpp ! distribued memory computing 31 USE lbclnk ! ocean lateral boundary condition (or mpp link) 32 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 31 USE lbclnk ! ocean lateral boundary condition (or mpp link) 32 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 33 33 34 34 IMPLICIT NONE … … 36 36 37 37 PUBLIC tra_adv_mus ! routine called by traadv.F90 38 38 39 39 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits 40 40 ! ! and in closed seas (orca 2 and 1 configurations) 41 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xind !: mixed upstream/centered index 42 42 43 43 LOGICAL :: l_trd ! flag to compute trends 44 44 LOGICAL :: l_ptr ! flag to compute poleward transport … … 50 50 !!---------------------------------------------------------------------- 51 51 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 52 !! $Id$ 52 !! $Id$ 53 53 !! Software governed by the CeCILL license (see ./LICENSE) 54 54 !!---------------------------------------------------------------------- … … 65 65 !! 66 66 !! ** Method : MUSCL scheme plus centered scheme at ocean boundaries 67 !! ld_msc_ups=T : 67 !! ld_msc_ups=T : 68 68 !! 69 69 !! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends … … 134 134 ! !-- first guess of the slopes 135 135 zwx(:,:,jpk) = 0._wp ! bottom values 136 zwy(:,:,jpk) = 0._wp 136 zwy(:,:,jpk) = 0._wp 137 137 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 138 138 zwx(ji,jj,jk) = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) … … 188 188 CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kbb) ) 189 189 END IF 190 ! ! "Poleward" heat and salt transports 190 ! ! "Poleward" heat and salt transports 191 191 IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) ) 192 192 ! ! heat transport -
NEMO/trunk/src/OCE/TRA/traadv_qck.F90
r13982 r14072 19 19 USE trc_oce ! share passive tracers/Ocean variables 20 20 USE trd_oce ! trends: ocean variables 21 USE trdtra ! trends manager: tracers 21 USE trdtra ! trends manager: tracers 22 22 USE diaptr ! poleward transport diagnostics 23 23 USE iom … … 26 26 USE lib_mpp ! distribued memory computing 27 27 USE lbclnk ! ocean lateral boundary condition (or mpp link) 28 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 28 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 29 29 30 30 IMPLICIT NONE … … 112 112 ! 113 113 ! ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme 114 CALL tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs ) 115 CALL tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs ) 114 CALL tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs ) 115 CALL tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs ) 116 116 117 117 ! ! vertical fluxes are computed with the 2nd order centered scheme … … 142 142 DO jn = 1, kjpt ! tracer loop 143 143 ! ! =========== 144 zfu(:,:,:) = 0._wp ; zfc(:,:,:) = 0._wp 145 zfd(:,:,:) = 0._wp ; zwx(:,:,:) = 0._wp 144 zfu(:,:,:) = 0._wp ; zfc(:,:,:) = 0._wp 145 zfd(:,:,:) = 0._wp ; zwx(:,:,:) = 0._wp 146 146 ! 147 147 !!gm why not using a SHIFT instruction... … … 151 151 END_3D 152 152 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 153 153 154 154 ! 155 155 ! Horizontal advective fluxes 156 156 ! --------------------------- 157 157 DO_3D( 0, 0, nn_hls-1, 0, 1, jpkm1 ) 158 zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 159 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk) ! FU in the x-direction for T 158 zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 159 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk) ! FU in the x-direction for T 160 160 END_3D 161 161 ! 162 162 DO_3D( 0, 0, nn_hls-1, 0, 1, jpkm1 ) 163 zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 163 zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 164 164 zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 165 165 zwx(ji,jj,jk) = ABS( pU(ji,jj,jk) ) * p2dt / zdx ! (0<zc_cfl<1 : Courant number on x-direction) … … 167 167 zfd(ji,jj,jk) = zdir * pt(ji+1,jj,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji ,jj,jk,jn,Kbb) ! FD in the x-direction for T 168 168 END_3D 169 !--- Lateral boundary conditions 169 !--- Lateral boundary conditions 170 170 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwx(:,:,:), 'T', 1.0_wp ) 171 171 … … 227 227 DO jn = 1, kjpt ! tracer loop 228 228 ! ! =========== 229 zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0 230 zfd(:,:,:) = 0.0 ; zwy(:,:,:) = 0.0 231 ! 232 DO jk = 1, jpkm1 233 ! 229 zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0 230 zfd(:,:,:) = 0.0 ; zwy(:,:,:) = 0.0 231 ! 232 DO jk = 1, jpkm1 233 ! 234 234 !--- Computation of the ustream and downstream value of the tracer and the mask 235 235 DO_2D( nn_hls-1, nn_hls-1, 0, 0 ) … … 241 241 END DO 242 242 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 243 243 244 244 ! 245 245 ! Horizontal advective fluxes … … 247 247 ! 248 248 DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 ) 249 zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 250 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk) ! FU in the x-direction for T 249 zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 250 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk) ! FU in the x-direction for T 251 251 END_3D 252 252 ! 253 253 DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 ) 254 zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 254 zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 255 255 zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) 256 256 zwy(ji,jj,jk) = ABS( pV(ji,jj,jk) ) * p2dt / zdx ! (0<zc_cfl<1 : Courant number on x-direction) … … 259 259 END_3D 260 260 261 !--- Lateral boundary conditions 261 !--- Lateral boundary conditions 262 262 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp ) 263 263 … … 328 328 IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean) 329 329 DO_2D( 0, 0, 0, 0 ) 330 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) ! linear free surface 330 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) ! linear free surface 331 331 END_2D 332 332 ELSE ! no ocean cavities (only ocean surface) … … 354 354 !! ** Purpose : Computation of advective flux with Quickest scheme 355 355 !! 356 !! ** Method : 356 !! ** Method : 357 357 !!---------------------------------------------------------------------- 358 358 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: pfu ! second upwind point … … 361 361 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) :: puc ! input as Courant number ; output as flux 362 362 !! 363 INTEGER :: ji, jj, jk ! dummy loop indices 364 REAL(wp) :: zcoef1, zcoef2, zcoef3 ! local scalars 363 INTEGER :: ji, jj, jk ! dummy loop indices 364 REAL(wp) :: zcoef1, zcoef2, zcoef3 ! local scalars 365 365 REAL(wp) :: zc, zcurv, zfho ! - - 366 366 !---------------------------------------------------------------------- … … 372 372 zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) ) 373 373 zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv 374 zfho = zcoef1 - zcoef2 - zcoef3 ! phi_f QUICKEST 374 zfho = zcoef1 - zcoef2 - zcoef3 ! phi_f QUICKEST 375 375 ! 376 376 zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk) … … 378 378 zcoef3 = ABS( zcurv ) 379 379 IF( zcoef3 >= zcoef2 ) THEN 380 zfho = pfc(ji,jj,jk) 380 zfho = pfc(ji,jj,jk) 381 381 ELSE 382 382 zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) ) ! phi_REF 383 383 IF( zcoef1 >= 0. ) THEN 384 zfho = MAX( pfc(ji,jj,jk), zfho ) 385 zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 384 zfho = MAX( pfc(ji,jj,jk), zfho ) 385 zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 386 386 ELSE 387 zfho = MIN( pfc(ji,jj,jk), zfho ) 388 zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 387 zfho = MIN( pfc(ji,jj,jk), zfho ) 388 zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 389 389 ENDIF 390 390 ENDIF -
NEMO/trunk/src/OCE/TRA/traadv_ubs.F90
r13982 r14072 10 10 !!---------------------------------------------------------------------- 11 11 !! tra_adv_ubs : update the tracer trend with the horizontal 12 !! advection trends using a third order biaised scheme 12 !! advection trends using a third order biaised scheme 13 13 !!---------------------------------------------------------------------- 14 14 USE oce ! ocean dynamics and active tracers … … 16 16 USE trc_oce ! share passive tracers/Ocean variables 17 17 USE trd_oce ! trends: ocean variables 18 USE traadv_fct ! acces to routine interp_4th_cpt 19 USE trdtra ! trends manager: tracers 18 USE traadv_fct ! acces to routine interp_4th_cpt 19 USE trdtra ! trends manager: tracers 20 20 USE diaptr ! poleward transport diagnostics 21 21 USE diaar5 ! AR5 diagnostics … … 25 25 USE lib_mpp ! massively parallel library 26 26 USE lbclnk ! ocean lateral boundary condition (or mpp link) 27 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 27 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 28 28 29 29 IMPLICIT NONE … … 51 51 !!---------------------------------------------------------------------- 52 52 !! *** ROUTINE tra_adv_ubs *** 53 !! 53 !! 54 54 !! ** Purpose : Compute the now trend due to the advection of tracers 55 55 !! and add it to the general trend of passive tracer equations. … … 60 60 !! For example the i-component of the advective fluxes are given by : 61 61 !! ! e2u e3u un ( mi(Tn) - zltu(i ) ) if un(i) >= 0 62 !! ztu = ! or 62 !! ztu = ! or 63 63 !! ! e2u e3u un ( mi(Tn) - zltu(i+1) ) if un(i) < 0 64 64 !! where zltu is the second derivative of the before temperature field: 65 65 !! zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ] 66 !! This results in a dissipatively dominant (i.e. hyper-diffusive) 67 !! truncation error. The overall performance of the advection scheme 68 !! is similar to that reported in (Farrow and Stevens, 1995). 66 !! This results in a dissipatively dominant (i.e. hyper-diffusive) 67 !! truncation error. The overall performance of the advection scheme 68 !! is similar to that reported in (Farrow and Stevens, 1995). 69 69 !! For stability reasons, the first term of the fluxes which corresponds 70 !! to a second order centered scheme is evaluated using the now velocity 71 !! (centered in time) while the second term which is the diffusive part 72 !! of the scheme, is evaluated using the before velocity (forward in time). 70 !! to a second order centered scheme is evaluated using the now velocity 71 !! (centered in time) while the second term which is the diffusive part 72 !! of the scheme, is evaluated using the before velocity (forward in time). 73 73 !! Note that UBS is not positive. Do not use it on passive tracers. 74 74 !! On the vertical, the advection is evaluated using a FCT scheme, 75 !! as the UBS have been found to be too diffusive. 76 !! kn_ubs_v argument controles whether the FCT is based on 77 !! a 2nd order centrered scheme (kn_ubs_v=2) or on a 4th order compact 75 !! as the UBS have been found to be too diffusive. 76 !! kn_ubs_v argument controles whether the FCT is based on 77 !! a 2nd order centrered scheme (kn_ubs_v=2) or on a 4th order compact 78 78 !! scheme (kn_ubs_v=4). 79 79 !! … … 82 82 !! - poleward advective heat and salt transport (ln_diaptr=T) 83 83 !! 84 !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404. 84 !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404. 85 85 !! Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731�1741. 86 86 !!---------------------------------------------------------------------- … … 125 125 DO jn = 1, kjpt ! tracer loop 126 126 ! ! =========== 127 ! 127 ! 128 128 DO jk = 1, jpkm1 !== horizontal laplacian of before tracer ==! 129 129 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) ! First derivative (masked gradient) … … 138 138 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zcoef 139 139 END_2D 140 ! 141 END DO 140 ! 141 END DO 142 142 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_ubs', zltu, 'T', 1.0_wp, zltv, 'T', 1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 143 ! 143 ! 144 144 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) !== Horizontal advective fluxes ==! (UBS) 145 145 zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) ) ! upstream transport (x2) … … 166 166 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 167 167 END_2D 168 ! 168 ! 169 169 END DO 170 170 ! … … 177 177 CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztv, pV, pt(:,:,:,jn,Kmm) ) 178 178 END IF 179 ! 179 ! 180 180 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 181 181 IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) ) … … 188 188 SELECT CASE( kn_ubs_v ) ! select the vertical advection scheme 189 189 ! 190 CASE( 2 ) ! 2nd order FCT 191 ! 190 CASE( 2 ) ! 2nd order FCT 191 ! 192 192 IF( l_trd ) THEN 193 193 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) … … 205 205 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface 206 206 DO_2D( 1, 1, 1, 1 ) 207 ztw(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface 207 ztw(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface 208 208 END_2D 209 209 ELSE ! no cavities: only at the ocean surface … … 217 217 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) & 218 218 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 219 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztak 219 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztak 220 220 zti(ji,jj,jk) = ( pt(ji,jj,jk,jn,Kbb) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 221 221 END_3D … … 266 266 !!--------------------------------------------------------------------- 267 267 !! *** ROUTINE nonosc_z *** 268 !! 269 !! ** Purpose : compute monotonic tracer fluxes from the upstream 270 !! scheme and the before field by a nonoscillatory algorithm 268 !! 269 !! ** Purpose : compute monotonic tracer fluxes from the upstream 270 !! scheme and the before field by a nonoscillatory algorithm 271 271 !! 272 272 !! ** Method : ... ??? -
NEMO/trunk/src/OCE/TRA/traatf.F90
r14045 r14072 26 26 !!---------------------------------------------------------------------- 27 27 USE oce ! ocean dynamics and tracers variables 28 USE dom_oce ! ocean space and time domain variables 28 USE dom_oce ! ocean space and time domain variables 29 29 USE sbc_oce ! surface boundary condition: ocean 30 30 USE sbcrnf ! river runoffs … … 33 33 USE domvvl ! variable volume 34 34 USE trd_oce ! trends: ocean variables 35 USE trdtra ! trends manager: tracers 35 USE trdtra ! trends manager: tracers 36 36 USE traqsr ! penetrative solar radiation (needed for nksr) 37 37 USE phycst ! physical constant … … 70 70 !! *** ROUTINE traatf *** 71 71 !! 72 !! ** Purpose : Apply the boundary condition on the after temperature 72 !! ** Purpose : Apply the boundary condition on the after temperature 73 73 !! and salinity fields and add the Asselin time filter on now fields. 74 !! 75 !! ** Method : At this stage of the computation, ta and sa are the 74 !! 75 !! ** Method : At this stage of the computation, ta and sa are the 76 76 !! after temperature and salinity as the time stepping has 77 77 !! been performed in trazdf_imp or trazdf_exp module. 78 78 !! 79 !! - Apply lateral boundary conditions on (ta,sa) 80 !! at the local domain boundaries through lbc_lnk call, 81 !! at the one-way open boundaries (ln_bdy=T), 79 !! - Apply lateral boundary conditions on (ta,sa) 80 !! at the local domain boundaries through lbc_lnk call, 81 !! at the one-way open boundaries (ln_bdy=T), 82 82 !! at the AGRIF zoom boundaries (lk_agrif=T) 83 83 !! … … 89 89 INTEGER , INTENT(in ) :: kt ! ocean time-step index 90 90 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level indices 91 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers 91 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers 92 92 !! 93 93 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 105 105 106 106 ! Update after tracer on domain lateral boundaries 107 ! 107 ! 108 108 #if defined key_agrif 109 109 CALL Agrif_tra ! AGRIF zoom boundaries … … 113 113 ! 114 114 IF( ln_bdy ) CALL bdy_tra( kt, Kbb, pts, Kaa ) ! BDY open boundaries 115 115 116 116 ! trends computation initialisation 117 IF( l_trdtra ) THEN 117 IF( l_trdtra ) THEN 118 118 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 119 119 ztrdt(:,:,:) = 0._wp 120 120 ztrds(:,:,:) = 0._wp 121 IF( ln_traldf_iso ) THEN ! diagnose the "pure" Kz diffusive trend 121 IF( ln_traldf_iso ) THEN ! diagnose the "pure" Kz diffusive trend 122 122 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_zdfp, ztrdt ) 123 123 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_zdfp, ztrds ) 124 124 ENDIF 125 ! total trend for the non-time-filtered variables. 125 ! total trend for the non-time-filtered variables. 126 126 zfact = 1.0 / rn_Dt 127 127 ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from pts(Kmm) terms … … 133 133 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_tot, ztrds ) 134 134 IF( ln_linssh ) THEN ! linear sea surface height only 135 ! Store now fields before applying the Asselin filter 135 ! Store now fields before applying the Asselin filter 136 136 ! in order to calculate Asselin filter trend later. 137 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kmm) 137 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kmm) 138 138 ztrds(:,:,:) = pts(:,:,:,jp_sal,Kmm) 139 139 ENDIF 140 140 ENDIF 141 141 142 IF( l_1st_euler ) THEN ! Euler time-stepping 142 IF( l_1st_euler ) THEN ! Euler time-stepping 143 143 ! 144 144 IF (l_trdtra .AND. .NOT. ln_linssh ) THEN ! Zero Asselin filter contribution must be explicitly written out since for vvl … … 152 152 ELSE ! Leap-Frog + Asselin filter time stepping 153 153 ! 154 IF( ln_linssh ) THEN ; CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000, 'TRA', pts, jpts ) ! linear free surface 154 IF( ln_linssh ) THEN ; CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000, 'TRA', pts, jpts ) ! linear free surface 155 155 ELSE ; CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nit000, rn_Dt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts ) ! non-linear free surface 156 156 ENDIF 157 157 ! 158 CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kmm) , 'T', 1.0_wp, pts(:,:,:,jp_sal,Kmm) , 'T', 1.0_wp ) 159 160 ENDIF 161 ! 162 IF( l_trdtra .AND. ln_linssh ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt 158 CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kmm) , 'T', 1.0_wp, pts(:,:,:,jp_sal,Kmm) , 'T', 1.0_wp ) 159 160 ENDIF 161 ! 162 IF( l_trdtra .AND. ln_linssh ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt 163 163 DO jk = 1, jpkm1 164 164 ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kmm) - ztrdt(:,:,jk) ) * r1_Dt … … 184 184 !! 185 185 !! ** Purpose : fixed volume: apply the Asselin time filter to the "now" field 186 !! 186 !! 187 187 !! ** Method : - Apply a Asselin time filter on now fields. 188 188 !! … … 209 209 ! 210 210 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 211 ztn = pt(ji,jj,jk,jn,Kmm) 211 ztn = pt(ji,jj,jk,jn,Kmm) 212 212 ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb) ! time laplacian on tracers 213 213 ! … … 224 224 !! *** ROUTINE tra_atf_vvl *** 225 225 !! 226 !! ** Purpose : Time varying volume: apply the Asselin time filter 227 !! 226 !! ** Purpose : Time varying volume: apply the Asselin time filter 227 !! 228 228 !! ** Method : - Apply a thickness weighted Asselin time filter on now fields. 229 229 !! pt(Kmm) = ( e3t_Kmm*pt(Kmm) + rn_atfp*[ e3t_Kbb*pt(Kbb) - 2 e3t_Kmm*pt(Kmm) + e3t_Kaa*pt(Kaa) ] ) … … 255 255 ENDIF 256 256 ! 257 IF( cdtype == 'TRA' ) THEN 257 IF( cdtype == 'TRA' ) THEN 258 258 ll_traqsr = ln_traqsr ! active tracers case and solar penetration 259 259 ll_rnf = ln_rnf ! active tracers case and river runoffs … … 261 261 ELSE ! passive tracers case 262 262 ll_traqsr = .FALSE. ! NO solar penetration 263 ll_rnf = .FALSE. ! NO river runoffs ???? !!gm BUG ? 264 ll_isf = .FALSE. ! NO ice shelf melting/freezing !!gm BUG ?? 263 ll_rnf = .FALSE. ! NO river runoffs ???? !!gm BUG ? 264 ll_isf = .FALSE. ! NO ice shelf melting/freezing !!gm BUG ?? 265 265 ENDIF 266 266 ! … … 272 272 zfact1 = rn_atfp * p2dt 273 273 zfact2 = zfact1 * r1_rho0 274 DO jn = 1, kjpt 274 DO jn = 1, kjpt 275 275 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 276 276 ze3t_b = e3t(ji,jj,jk,Kbb) … … 289 289 ! 290 290 ! Add asselin correction on scale factors: 291 zscale = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) 292 ze3t_f = ze3t_f - zfact2 * zscale * ( emp_b(ji,jj) - emp(ji,jj) ) 293 IF ( ll_rnf ) ze3t_f = ze3t_f + zfact2 * zscale * ( rnf_b(ji,jj) - rnf(ji,jj) ) 291 zscale = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) 292 ze3t_f = ze3t_f - zfact2 * zscale * ( emp_b(ji,jj) - emp(ji,jj) ) 293 IF ( ll_rnf ) ze3t_f = ze3t_f + zfact2 * zscale * ( rnf_b(ji,jj) - rnf(ji,jj) ) 294 294 IF ( ll_isf ) THEN 295 295 IF ( ln_isfcav_mlt ) ze3t_f = ze3t_f - zfact2 * zscale * ( fwfisf_cav_b(ji,jj) - fwfisf_cav(ji,jj) ) … … 297 297 ENDIF 298 298 ! 299 IF( jk == mikt(ji,jj) ) THEN ! first level 299 IF( jk == mikt(ji,jj) ) THEN ! first level 300 300 ztc_f = ztc_f - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) ) 301 301 ENDIF 302 302 ! 303 303 ! solar penetration (temperature only) 304 IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr ) & 305 & ztc_f = ztc_f - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 304 IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr ) & 305 & ztc_f = ztc_f - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 306 306 ! 307 307 ! 308 308 IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) ) & 309 & ztc_f = ztc_f - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 309 & ztc_f = ztc_f - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 310 310 & * e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj) 311 311 … … 321 321 & * e3t(ji,jj,jk,Kmm) / rhisf_tbl_cav(ji,jj) 322 322 END IF 323 ! level partially include in Losch_2008 ice shelf boundary layer 323 ! level partially include in Losch_2008 ice shelf boundary layer 324 324 IF ( jk == misfkb_cav(ji,jj) ) THEN 325 325 ztc_f = ztc_f - zfact1 * ( risf_cav_tsc(ji,jj,jn) - risf_cav_tsc_b(ji,jj,jn) ) & … … 335 335 & * e3t(ji,jj,jk,Kmm) / rhisf_tbl_par(ji,jj) 336 336 END IF 337 ! level partially include in Losch_2008 ice shelf boundary layer 337 ! level partially include in Losch_2008 ice shelf boundary layer 338 338 IF ( jk == misfkb_par(ji,jj) ) THEN 339 339 ztc_f = ztc_f - zfact1 * ( risf_par_tsc(ji,jj,jn) - risf_par_tsc_b(ji,jj,jn) ) & … … 364 364 ! 365 365 END_3D 366 ! 366 ! 367 367 END DO 368 368 ! 369 369 IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) ) THEN 370 IF( l_trdtra .AND. cdtype == 'TRA' ) THEN 370 IF( l_trdtra .AND. cdtype == 'TRA' ) THEN 371 371 CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) ) 372 372 CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) ) -
NEMO/trunk/src/OCE/TRA/traatf_qco.F90
r14053 r14072 100 100 IF(lwp) WRITE(numout,*) '~~~~~~~' 101 101 ENDIF 102 !!st Update after tracer on domain lateral boundaries as been removed outside 102 !!st Update after tracer on domain lateral boundaries as been removed outside 103 103 104 104 ! trends computation initialisation -
NEMO/trunk/src/OCE/TRA/trabbc.F90
r13982 r14072 12 12 13 13 !!---------------------------------------------------------------------- 14 !! tra_bbc : update the tracer trend at ocean bottom 14 !! tra_bbc : update the tracer trend at ocean bottom 15 15 !! tra_bbc_init : initialization of geothermal heat flux trend 16 16 !!---------------------------------------------------------------------- … … 19 19 USE phycst ! physical constants 20 20 USE trd_oce ! trends: ocean variables 21 USE trdtra ! trends manager: tracers 21 USE trdtra ! trends manager: tracers 22 22 ! 23 23 USE in_out_manager ! I/O manager 24 USE iom ! xIOS 24 USE iom ! xIOS 25 25 USE fldread ! read input fields 26 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 43 43 44 44 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_qgh ! structure of input qgh (file informations, fields read) 45 45 46 46 !! * Substitutions 47 47 # include "do_loop_substitute.h90" … … 58 58 !! *** ROUTINE tra_bbc *** 59 59 !! 60 !! ** Purpose : Compute the bottom boundary contition on temperature 61 !! associated with geothermal heating and add it to the 60 !! ** Purpose : Compute the bottom boundary contition on temperature 61 !! associated with geothermal heating and add it to the 62 62 !! general trend of temperature equations. 63 63 !! 64 !! ** Method : The geothermal heat flux set to its constant value of 64 !! ** Method : The geothermal heat flux set to its constant value of 65 65 !! 86.4 mW/m2 (Stein and Stein 1992, Huang 1999). 66 66 !! The temperature trend associated to this heat flux through the … … 135 135 CHARACTER(len=256) :: cn_dir ! Root directory for location of ssr files 136 136 !! 137 NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir 137 NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir 138 138 !!---------------------------------------------------------------------- 139 139 ! -
NEMO/trunk/src/OCE/TRA/trabbl.F90
r13982 r14072 31 31 USE trdtra ! trends: active tracers 32 32 ! 33 USE iom ! IOM library 33 USE iom ! IOM library 34 34 USE in_out_manager ! I/O manager 35 35 USE lbclnk ! ocean lateral boundary conditions 36 36 USE prtctl ! Print control 37 37 USE timing ! Timing 38 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 38 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 39 39 40 40 IMPLICIT NONE … … 200 200 zptb(ji,jj) = pt(ji,jj,ik,jn) ! bottom before T and S 201 201 END_2D 202 ! 202 ! 203 203 DO_2D( 0, 0, 0, 0 ) ! Compute the trend 204 204 ik = mbkt(ji,jj) ! bottom T-level index … … 399 399 za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at u-point 400 400 zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 401 ! ! 2*masked bottom density gradient 401 ! ! 2*masked bottom density gradient 402 402 zgdrho = ( za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) ) & 403 403 - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) ) ) * umask(ji,jj,1) … … 523 523 END_2D 524 524 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 525 zmbku(:,:) = REAL( mbku_d(:,:), wp ) ; zmbkv(:,:) = REAL( mbkv_d(:,:), wp ) 526 CALL lbc_lnk_multi( 'trabbl', zmbku,'U',1.0_wp, zmbkv,'V',1.0_wp) 525 zmbku(:,:) = REAL( mbku_d(:,:), wp ) ; zmbkv(:,:) = REAL( mbkv_d(:,:), wp ) 526 CALL lbc_lnk_multi( 'trabbl', zmbku,'U',1.0_wp, zmbkv,'V',1.0_wp) 527 527 mbku_d(:,:) = MAX( INT( zmbku(:,:) ), 1 ) ; mbkv_d(:,:) = MAX( NINT( zmbkv(:,:) ), 1 ) 528 528 ! -
NEMO/trunk/src/OCE/TRA/tradmp.F90
r13982 r14072 11 11 !! NEMO 1.0 ! 2002-08 (G. Madec, E. Durand) free form + modules 12 12 !! 3.2 ! 2009-08 (G. Madec, C. Talandier) DOCTOR norm for namelist parameter 13 !! 3.3 ! 2010-06 (C. Ethe, G. Madec) merge TRA-TRC 13 !! 3.3 ! 2010-06 (C. Ethe, G. Madec) merge TRA-TRC 14 14 !! 3.4 ! 2011-04 (G. Madec, C. Ethe) Merge of dtatem and dtasal + suppression of CPP keys 15 15 !! 3.6 ! 2015-06 (T. Graham) read restoring coefficient in a file … … 26 26 USE c1d ! 1D vertical configuration 27 27 USE trd_oce ! trends: ocean variables 28 USE trdtra ! trends manager: tracers 28 USE trdtra ! trends manager: tracers 29 29 USE zdf_oce ! ocean: vertical physics 30 30 USE phycst ! physical constants … … 55 55 !!---------------------------------------------------------------------- 56 56 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 57 !! $Id$ 57 !! $Id$ 58 58 !! Software governed by the CeCILL license (see ./LICENSE) 59 59 !!---------------------------------------------------------------------- … … 75 75 !!---------------------------------------------------------------------- 76 76 !! *** ROUTINE tra_dmp *** 77 !! 77 !! 78 78 !! ** Purpose : Compute the tracer trend due to a newtonian damping 79 79 !! of the tracer field towards given data field and add it to the 80 80 !! general tracer trends. 81 81 !! 82 !! ** Method : Newtonian damping towards t_dta and s_dta computed 82 !! ** Method : Newtonian damping towards t_dta and s_dta computed 83 83 !! and add to the general tracer trends: 84 84 !! ta = ta + resto * (t_dta - tb) … … 158 158 !!---------------------------------------------------------------------- 159 159 !! *** ROUTINE tra_dmp_init *** 160 !! 161 !! ** Purpose : Initialization for the newtonian damping 160 !! 161 !! ** Purpose : Initialization for the newtonian damping 162 162 !! 163 163 !! ** Method : read the namtra_dmp namelist and check the parameters 164 164 !!---------------------------------------------------------------------- 165 INTEGER :: ios, imask ! local integers 165 INTEGER :: ios, imask ! local integers 166 166 ! 167 167 NAMELIST/namtra_dmp/ ln_tradmp, nn_zdmp, cn_resto -
NEMO/trunk/src/OCE/TRA/traisf.F90
r13982 r14072 35 35 !!---------------------------------------------------------------------- 36 36 !! *** ROUTINE tra_isf *** 37 !! 37 !! 38 38 !! ** Purpose : Compute the temperature trend due to the ice shelf melting (qhoce + qhc) 39 39 !! … … 65 65 ! 66 66 ! Dynamical stability at start up after change in under ice shelf cavity geometry is achieve by correcting the divergence. 67 ! This is achieved by applying a volume flux in order to keep the horizontal divergence after remapping 67 ! This is achieved by applying a volume flux in order to keep the horizontal divergence after remapping 68 68 ! the same as at the end of the latest time step. So correction need to be apply at nit000 (euler time step) and 69 69 ! half of it at nit000+1 (leap frog time step). … … 95 95 !! *** Purpose : Compute the temperature trend due to the ice shelf melting (qhoce + qhc) for cav or par case 96 96 !! 97 !! *** Action :: Update pts(:,:,:,:,Krhs) with the surface boundary condition trend 97 !! *** Action :: Update pts(:,:,:,:,Krhs) with the surface boundary condition trend 98 98 !! 99 99 !!---------------------------------------------------------------------- … … 104 104 REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: ptsc , ptsc_b 105 105 !!---------------------------------------------------------------------- 106 INTEGER :: ji,jj,jk ! loop index 106 INTEGER :: ji,jj,jk ! loop index 107 107 INTEGER :: ikt, ikb ! top and bottom level of the tbl 108 108 REAL(wp), DIMENSION(A2D(nn_hls)) :: ztc ! total ice shelf tracer trend … … 125 125 END DO 126 126 ! 127 ! level partially include in ice shelf boundary layer 127 ! level partially include in ice shelf boundary layer 128 128 pts(ji,jj,ikb,jp_tem) = pts(ji,jj,ikb,jp_tem) + ztc(ji,jj) * pfrac(ji,jj) 129 129 ! … … 136 136 !! *** ROUTINE tra_isf_cpl *** 137 137 !! 138 !! *** Action :: Update pts(:,:,:,:,Krhs) with the ice shelf coupling trend 138 !! *** Action :: Update pts(:,:,:,:,Krhs) with the ice shelf coupling trend 139 139 !! 140 140 !!---------------------------------------------------------------------- -
NEMO/trunk/src/OCE/TRA/traldf.F90
r13982 r14072 2 2 !!====================================================================== 3 3 !! *** MODULE traldf *** 4 !! Ocean Active tracers : lateral diffusive trends 4 !! Ocean Active tracers : lateral diffusive trends 5 5 !!===================================================================== 6 6 !! History : 9.0 ! 2005-11 (G. Madec) Original code 7 !! NEMO 3.0 ! 2008-01 (C. Ethe, G. Madec) merge TRC-TRA 7 !! NEMO 3.0 ! 2008-01 (C. Ethe, G. Madec) merge TRC-TRA 8 8 !! 3.7 ! 2013-12 (G. Madec) remove the optional computation from T & S anomaly profiles and traldf_bilapg 9 9 !! - ! 2013-12 (F. Lemarie, G. Madec) triad operator (Griffies) + Method of Stabilizing Correction … … 37 37 PRIVATE 38 38 39 PUBLIC tra_ldf ! called by step.F90 40 PUBLIC tra_ldf_init ! called by nemogcm.F90 39 PUBLIC tra_ldf ! called by step.F90 40 PUBLIC tra_ldf_init ! called by nemogcm.F90 41 41 42 42 !!---------------------------------------------------------------------- 43 43 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 44 !! $Id$ 44 !! $Id$ 45 45 !! Software governed by the CeCILL license (see ./LICENSE) 46 46 !!---------------------------------------------------------------------- … … 50 50 !!---------------------------------------------------------------------- 51 51 !! *** ROUTINE tra_ldf *** 52 !! 52 !! 53 53 !! ** Purpose : compute the lateral ocean tracer physics. 54 54 !!---------------------------------------------------------------------- … … 120 120 !!---------------------------------------------------------------------- 121 121 !! *** ROUTINE tra_ldf_init *** 122 !! 122 !! 123 123 !! ** Purpose : Choice of the operator for the lateral tracer diffusion 124 124 !! 125 125 !! ** Method : set nldf_tra from the namtra_ldf logicals 126 126 !!---------------------------------------------------------------------- 127 INTEGER :: ioptio, ierr ! temporary integers 127 INTEGER :: ioptio, ierr ! temporary integers 128 128 !!---------------------------------------------------------------------- 129 129 ! -
NEMO/trunk/src/OCE/TRA/traldf_iso.F90
r13982 r14072 15 15 !!---------------------------------------------------------------------- 16 16 !! tra_ldf_iso : update the tracer trend with the horizontal component of a iso-neutral laplacian operator 17 !! and with the vertical part of the isopycnal or geopotential s-coord. operator 17 !! and with the vertical part of the isopycnal or geopotential s-coord. operator 18 18 !!---------------------------------------------------------------------- 19 19 USE oce ! ocean dynamics and active tracers … … 79 79 !! *** ROUTINE tra_ldf_iso *** 80 80 !! 81 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive 82 !! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 81 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive 82 !! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 83 83 !! add it to the general trend of tracer equation. 84 84 !! 85 !! ** Method : The horizontal component of the lateral diffusive trends 85 !! ** Method : The horizontal component of the lateral diffusive trends 86 86 !! is provided by a 2nd order operator rotated along neural or geopo- 87 87 !! tential surfaces to which an eddy induced advection can be added … … 94 94 !! 95 95 !! 2nd part : horizontal fluxes of the lateral mixing operator 96 !! ======== 96 !! ======== 97 97 !! zftu = pahu e2u*e3u/e1u di[ tb ] 98 98 !! - pahu e2u*uslp dk[ mi(mk(tb)) ] … … 165 165 ELSE ; zsign = -1._wp 166 166 ENDIF 167 167 168 168 !!---------------------------------------------------------------------- 169 169 !! 0 - calculate ah_wslp2 and akz … … 223 223 DO jn = 1, kjpt ! tracer loop 224 224 ! ! =========== 225 ! 226 !!---------------------------------------------------------------------- 227 !! I - masked horizontal derivative 225 ! 226 !!---------------------------------------------------------------------- 227 !! I - masked horizontal derivative 228 228 !!---------------------------------------------------------------------- 229 229 !!gm : bug.... why (x,:,:)? (1,jpj,:) and (jpi,1,:) should be sufficient.... … … 232 232 !!end 233 233 234 ! Horizontal tracer gradient 234 ! Horizontal tracer gradient 235 235 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 236 236 zdit(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) … … 239 239 IF( ln_zps ) THEN ! botton and surface ocean correction of the horizontal gradient 240 240 DO_2D( 1, 0, 1, 0 ) ! bottom correction (partial bottom cell) 241 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 241 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 242 242 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 243 243 END_2D 244 244 IF( ln_isfcav ) THEN ! first wet level beneath a cavity 245 245 DO_2D( 1, 0, 1, 0 ) 246 IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 247 IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 246 IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 247 IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 248 248 END_2D 249 249 ENDIF … … 283 283 zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & 284 284 & + zcof2 * ( zdkt (ji,jj+1) + zdk1t(ji,jj) & 285 & + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk) 285 & + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk) 286 286 END_2D 287 287 ! … … 291 291 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 292 292 END_2D 293 END DO ! End of slab 293 END DO ! End of slab 294 294 295 295 !!---------------------------------------------------------------------- … … 301 301 ! ! Surface and bottom vertical fluxes set to zero 302 302 ztfw(:,:, 1 ) = 0._wp ; ztfw(:,:,jpk) = 0._wp 303 303 304 304 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! interior (2=<jk=<jpk-1) 305 305 ! … … 330 330 END_3D 331 331 ! 332 ELSE ! bilaplacian 332 ELSE ! bilaplacian 333 333 SELECT CASE( kpass ) 334 334 CASE( 1 ) ! 1st pass : eddy coef = ah_wslp2 … … 346 346 END SELECT 347 347 ENDIF 348 ! 348 ! 349 349 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== Divergence of vertical fluxes added to pta ==! 350 350 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) & -
NEMO/trunk/src/OCE/TRA/traldf_lap_blp.F90
r13982 r14072 4 4 !! Ocean tracers: lateral diffusivity trend (laplacian and bilaplacian) 5 5 !!============================================================================== 6 !! History : 3.7 ! 2014-01 (G. Madec, S. Masson) Original code, re-entrant laplacian 6 !! History : 3.7 ! 2014-01 (G. Madec, S. Masson) Original code, re-entrant laplacian 7 7 !!---------------------------------------------------------------------- 8 8 … … 74 74 !!---------------------------------------------------------------------- 75 75 !! *** ROUTINE tra_ldf_lap *** 76 !! 77 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive 76 !! 77 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive 78 78 !! trend and add it to the general trend of tracer equation. 79 79 !! 80 80 !! ** Method : Second order diffusive operator evaluated using before 81 !! fields (forward time scheme). The horizontal diffusive trends of 81 !! fields (forward time scheme). The horizontal diffusive trends of 82 82 !! the tracer is given by: 83 83 !! difft = 1/(e1e2t*e3t) { di-1[ pahu e2u*e3u/e1u di(tb) ] … … 86 86 !! pt_rhs = pt_rhs + difft 87 87 !! 88 !! ** Action : - Update pt_rhs arrays with the before iso-level 88 !! ** Action : - Update pt_rhs arrays with the before iso-level 89 89 !! harmonic mixing trend. 90 90 !!---------------------------------------------------------------------- … … 139 139 ! ! =========== ! 140 140 DO jn = 1, kjpt ! tracer loop ! 141 ! ! =========== ! 142 ! 141 ! ! =========== ! 142 ! 143 143 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) !== First derivative (gradient) ==! 144 144 ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt(ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) ) … … 152 152 IF( ln_isfcav ) THEN ! top in ocean cavities only 153 153 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 154 IF( miku(ji,jj) > 1 ) ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn) 155 IF( mikv(ji,jj) > 1 ) ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn) 154 IF( miku(ji,jj) > 1 ) ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn) 155 IF( mikv(ji,jj) > 1 ) ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn) 156 156 END_2D 157 157 ENDIF … … 177 177 ! 178 178 END SUBROUTINE tra_ldf_lap_t 179 179 180 180 181 181 SUBROUTINE tra_ldf_blp( kt, Kmm, kit000, cdtype, pahu, pahv , & … … 184 184 !!---------------------------------------------------------------------- 185 185 !! *** ROUTINE tra_ldf_blp *** 186 !! 187 !! ** Purpose : Compute the before lateral tracer diffusive 186 !! 187 !! ** Purpose : Compute the before lateral tracer diffusive 188 188 !! trend and add it to the general trend of tracer equation. 189 189 !! … … 238 238 ! NOTE: [tiling-comms-merge] Needed for both nn_hls as tra_ldf_iso and tra_ldf_triad have not yet been adjusted to work with nn_hls = 2. In the zps case the lbc_lnk in zps_hde handles this, but in the zco case zlap always needs this lbc_lnk. I did try adjusting the bounds in tra_ldf_iso and tra_ldf_triad so this lbc_lnk was only needed for nn_hls = 1, but this was not correct and I did not have time to figure out why 239 239 CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1.0_wp ) ! Lateral boundary conditions (unchanged sign) 240 ! ! Partial top/bottom cell: GRADh( zlap ) 240 ! ! Partial top/bottom cell: GRADh( zlap ) 241 241 IF( ln_isfcav .AND. ln_zps ) THEN ; CALL zps_hde_isf( kt, Kmm, kjpt, zlap, zglu, zglv, zgui, zgvi ) ! both top & bottom 242 ELSEIF( ln_zps ) THEN ; CALL zps_hde ( kt, Kmm, kjpt, zlap, zglu, zglv ) ! only bottom 242 ELSEIF( ln_zps ) THEN ; CALL zps_hde ( kt, Kmm, kjpt, zlap, zglu, zglv ) ! only bottom 243 243 ENDIF 244 244 ! -
NEMO/trunk/src/OCE/TRA/traldf_triad.F90
r13982 r14072 145 145 ELSE ; zsign = -1._wp 146 146 ENDIF 147 ! 147 ! 148 148 !!---------------------------------------------------------------------- 149 149 !! 0 - calculate ah_wslp2, akz, and optionally zpsi_uw, zpsi_vw … … 175 175 END DO 176 176 ! 177 DO jp = 0, 1 ! j-k triads 177 DO jp = 0, 1 ! j-k triads 178 178 DO kp = 0, 1 179 179 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) … … 264 264 IF( ln_isfcav ) THEN ! top level (ocean cavities only) 265 265 DO_2D( 1, 0, 1, 0 ) 266 IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn) 267 IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn) 266 IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn) 267 IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn) 268 268 END_2D 269 269 ENDIF … … 392 392 & * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 393 393 END_3D 394 ELSE ! bilaplacian 394 ELSE ! bilaplacian 395 395 SELECT CASE( kpass ) 396 396 CASE( 1 ) ! 1st pass : eddy coef = ah_wslp2 … … 405 405 & + akz (ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) ) ) 406 406 END_3D 407 END SELECT 407 END SELECT 408 408 ENDIF 409 409 ! -
NEMO/trunk/src/OCE/TRA/tranpc.F90
r13982 r14072 97 97 IF( l_LB_debug ) THEN 98 98 ! Location of 1 known convection site to follow what's happening in the water column 99 ilc1 = 45 ; jlc1 = 3 ; ! ORCA2 4x4, Antarctic coast, more than 2 unstable portions in the water column... 99 ilc1 = 45 ; jlc1 = 3 ; ! ORCA2 4x4, Antarctic coast, more than 2 unstable portions in the water column... 100 100 nncpu = 1 ; ! the CPU domain contains the convection spot 101 klc1 = mbkt(ilc1,jlc1) ! bottom of the ocean for debug point... 101 klc1 = mbkt(ilc1,jlc1) ! bottom of the ocean for debug point... 102 102 ENDIF 103 103 ! … … 116 116 ! 117 117 IF( tmask(ji,jj,2) == 1 ) THEN ! At least 2 ocean points 118 ! ! consider one ocean column 118 ! ! consider one ocean column 119 119 zvts(:,jp_tem) = pts(ji,jj,:,jp_tem,Kaa) ! temperature 120 120 zvts(:,jp_sal) = pts(ji,jj,:,jp_sal,Kaa) ! salinity 121 121 ! 122 zvab(:,jp_tem) = zab(ji,jj,:,jp_tem) ! Alpha 123 zvab(:,jp_sal) = zab(ji,jj,:,jp_sal) ! Beta 124 zvn2(:) = zn2(ji,jj,:) ! N^2 122 zvab(:,jp_tem) = zab(ji,jj,:,jp_tem) ! Alpha 123 zvab(:,jp_sal) = zab(ji,jj,:,jp_sal) ! Beta 124 zvn2(:) = zn2(ji,jj,:) ! N^2 125 125 ! 126 126 IF( l_LB_debug ) THEN !LB debug: … … 128 128 IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE. 129 129 ! writing only if on CPU domain where conv region is: 130 lp_monitor_point = (narea == nncpu).AND.lp_monitor_point 130 lp_monitor_point = (narea == nncpu).AND.lp_monitor_point 131 131 ENDIF !LB debug end 132 132 ! … … 140 140 ! 141 141 jiter = jiter + 1 142 ! 142 ! 143 143 IF( jiter >= 400 ) EXIT 144 144 ! … … 155 155 ilayer = ilayer + 1 ! yet another instable portion of the water column found.... 156 156 ! 157 IF( lp_monitor_point ) THEN 157 IF( lp_monitor_point ) THEN 158 158 WRITE(numout,*) 159 159 IF( ilayer == 1 .AND. jiter == 1 ) THEN ! first time a column is spoted with an instability … … 195 195 zsum_beta = 0._wp 196 196 zsum_z = 0._wp 197 197 198 198 DO jk = ikup, ikbot ! Inside the instable (and overlying neutral) portion of the column 199 199 ! … … 204 204 zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz 205 205 zsum_z = zsum_z + zdz 206 ! 206 ! 207 207 IF( jk == ikbot ) EXIT ! avoid array-index overshoot in case ikbot = jpk, cause we're calling jk+1 next line 208 208 !! EXIT when we have reached the last layer that is instable (N2<0) or neutral (N2=0): 209 209 IF( zvn2(jk+1) > zn2_zero ) EXIT 210 210 END DO 211 211 212 212 ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative or neutral N2 213 213 IF( ikup == ikdown ) CALL ctl_stop( 'tra_npc : PROBLEM #2') … … 235 235 zvab(jk,jp_sal) = zbeta 236 236 END DO 237 238 237 238 239 239 !! Updating N2 in the relvant portion of the water column 240 240 !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion 241 241 !! => Need to re-compute N2! will use Alpha and Beta! 242 242 243 243 ikup = MAX(2,ikup) ! ikup can never be 1 ! 244 244 ik_low = MIN(ikdown+1,ikbot) ! we must go 1 point deeper than ikdown! 245 245 246 246 DO jk = ikup, ik_low ! we must go 1 point deeper than ikdown! 247 247 … … 263 263 264 264 END DO 265 265 266 266 ikp = MIN(ikdown+1,ikbot) 267 267 268 268 269 269 ENDIF !IF( zvn2(ikp) < 0. ) … … 275 275 276 276 IF( ikp /= ikbot ) CALL ctl_stop( 'tra_npc : PROBLEM #3') 277 277 278 278 ! ******* At this stage ikp == ikbot ! ******* 279 279 280 280 IF( ilayer > 0 ) THEN !! least an unstable layer has been found 281 281 ! -
NEMO/trunk/src/OCE/TRA/traqsr.F90
r14053 r14072 9 9 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 10 10 !! - ! 2005-11 (G. Madec) zco, zps, sco coordinate 11 !! 3.2 ! 2009-04 (G. Madec & NEMO team) 12 !! 3.6 ! 2012-05 (C. Rousset) store attenuation coef for use in ice model 11 !! 3.2 ! 2009-04 (G. Madec & NEMO team) 12 !! 3.6 ! 2012-05 (C. Rousset) store attenuation coef for use in ice model 13 13 !! 3.6 ! 2015-12 (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll 14 !! 3.7 ! 2015-11 (G. Madec, A. Coward) remove optimisation for fix volume 14 !! 3.7 ! 2015-11 (G. Madec, A. Coward) remove optimisation for fix volume 15 15 !!---------------------------------------------------------------------- 16 16 17 17 !!---------------------------------------------------------------------- 18 !! tra_qsr : temperature trend due to the penetration of solar radiation 19 !! tra_qsr_init : initialization of the qsr penetration 18 !! tra_qsr : temperature trend due to the penetration of solar radiation 19 !! tra_qsr_init : initialization of the qsr penetration 20 20 !!---------------------------------------------------------------------- 21 21 USE oce ! ocean dynamics and active tracers … … 45 45 ! !!* Namelist namtra_qsr: penetrative solar radiation 46 46 LOGICAL , PUBLIC :: ln_traqsr !: light absorption (qsr) flag 47 LOGICAL , PUBLIC :: ln_qsr_rgb !: Red-Green-Blue light absorption flag 47 LOGICAL , PUBLIC :: ln_qsr_rgb !: Red-Green-Blue light absorption flag 48 48 LOGICAL , PUBLIC :: ln_qsr_2bd !: 2 band light absorption flag 49 49 LOGICAL , PUBLIC :: ln_qsr_bio !: bio-model light absorption flag … … 54 54 ! 55 55 INTEGER , PUBLIC :: nksr !: levels below which the light cannot penetrate (depth larger than 391 m) 56 56 57 57 INTEGER, PARAMETER :: np_RGB = 1 ! R-G-B light penetration with constant Chlorophyll 58 58 INTEGER, PARAMETER :: np_RGBc = 2 ! R-G-B light penetration with Chlorophyll data … … 88 88 !! Considering the 2 wavebands case: 89 89 !! I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 90 !! The temperature trend associated with the solar radiation penetration 90 !! The temperature trend associated with the solar radiation penetration 91 91 !! is given by : zta = 1/e3t dk[ I ] / (rho0*Cp) 92 92 !! At the bottom, boudary condition for the radiation is no flux : 93 93 !! all heat which has not been absorbed in the above levels is put 94 94 !! in the last ocean level. 95 !! The computation is only done down to the level where 96 !! I(k) < 1.e-15 W/m2 (i.e. over the top nksr levels) . 95 !! The computation is only done down to the level where 96 !! I(k) < 1.e-15 W/m2 (i.e. over the top nksr levels) . 97 97 !! 98 98 !! ** Action : - update ta with the penetrative solar radiation trend … … 193 193 DO_2D( isj, iej, isi, iei ) 194 194 ! zlogc = log(zchl) 195 zlogc = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) ) 195 zlogc = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) ) 196 196 ! zc1 : log(zCze) = log (1.12 * zchl**0.803) 197 zc1 = 0.113328685307 + 0.803 * zlogc 197 zc1 = 0.113328685307 + 0.803 * zlogc 198 198 ! zc2 : log(zCtot) = log(40.6 * zchl**0.459) 199 zc2 = 3.703768066608 + 0.459 * zlogc 199 zc2 = 3.703768066608 + 0.459 * zlogc 200 200 ! zc3 : log(zze) = log(568.2 * zCtot**(-0.746)) 201 zc3 = 6.34247346942 - 0.746 * zc2 201 zc3 = 6.34247346942 - 0.746 * zc2 202 202 ! IF( log(zze) > log(102.) ) log(zze) = log(200.0 * zCtot**(-0.293)) 203 IF( zc3 > 4.62497281328 ) zc3 = 5.298317366548 - 0.293 * zc2 204 ! 203 IF( zc3 > 4.62497281328 ) zc3 = 5.298317366548 - 0.293 * zc2 204 ! 205 205 ze0(ji,jj) = zlogc ! ze0 = log(zchl) 206 206 ze1(ji,jj) = EXP( zc1 ) ! ze1 = zCze … … 208 208 ze3(ji,jj) = EXP( - zc3 ) ! ze3 = 1/zze 209 209 END_2D 210 210 211 211 ! 212 212 DO_3D( isj, iej, isi, iei, 1, nksr + 1 ) … … 230 230 ELSE !* constant chlorophyll 231 231 zchl = 0.05 232 ! NB. make sure constant value is such that: 232 ! NB. make sure constant value is such that: 233 233 zchl = MIN( 10. , MAX( 0.03, zchl ) ) 234 234 ! Convert chlorophyll value to attenuation coefficient look-up table index … … 245 245 ze2(ji,jj) = zcoef * qsr(ji,jj) 246 246 ze3(ji,jj) = zcoef * qsr(ji,jj) 247 ! store the surface SW radiation; re-use the surface ztmp3d array 247 ! store the surface SW radiation; re-use the surface ztmp3d array 248 248 ! since the surface attenuation coefficient is not used 249 249 ztmp3d(ji,jj,1) = qsr(ji,jj) … … 269 269 END_3D 270 270 ! 271 DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d ) 271 DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d ) 272 272 ! 273 273 CASE( np_2BD ) !== 2-bands fluxes ==! … … 278 278 zc0 = zz0 * EXP( -gdepw(ji,jj,jk ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk ,Kmm)*xsi1r ) 279 279 zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r ) 280 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 280 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 281 281 END_3D 282 282 ! … … 341 341 !! from two length scale of penetration (rn_si0,rn_si1) and a ratio 342 342 !! (rn_abs). These parameters are read in the namtra_qsr namelist. The 343 !! default values correspond to clear water (type I in Jerlov' 343 !! default values correspond to clear water (type I in Jerlov' 344 344 !! (1968) classification. 345 345 !! called by tra_qsr at the first timestep (nit000) … … 391 391 & ' 2 bands, 3 RGB bands or bio-model light penetration' ) 392 392 ! 393 IF( ln_qsr_rgb .AND. nn_chldta == 0 ) nqsr = np_RGB 393 IF( ln_qsr_rgb .AND. nn_chldta == 0 ) nqsr = np_RGB 394 394 IF( ln_qsr_rgb .AND. nn_chldta == 1 ) nqsr = np_RGBc 395 395 IF( ln_qsr_2bd ) nqsr = np_2BD … … 401 401 ! 402 402 SELECT CASE( nqsr ) 403 ! 403 ! 404 404 CASE( np_RGB , np_RGBc ) !== Red-Green-Blue light penetration ==! 405 ! 405 ! 406 406 IF(lwp) WRITE(numout,*) ' ==>>> R-G-B light penetration ' 407 407 ! 408 408 CALL trc_oce_rgb( rkrgb ) ! tabulated attenuation coef. 409 ! 409 ! 410 410 nksr = trc_oce_ext_lev( r_si2, 33._wp ) ! level of light extinction 411 411 ! … … 441 441 ! 442 442 CALL trc_oce_rgb( rkrgb ) ! tabulated attenuation coef. 443 ! 443 ! 444 444 nksr = trc_oce_ext_lev( r_si2, 33._wp ) ! level of light extinction 445 445 ! -
NEMO/trunk/src/OCE/TRA/trasbc.F90
r14053 r14072 9 9 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps 10 10 !! - ! 2010-09 (C. Ethe, G. Madec) Merge TRA-TRC 11 !! 3.6 ! 2014-11 (P. Mathiot) isf melting forcing 11 !! 3.6 ! 2014-11 (P. Mathiot) isf melting forcing 12 12 !! 4.1 ! 2019-09 (P. Mathiot) isf moved in traisf 13 13 !!---------------------------------------------------------------------- … … 21 21 USE phycst ! physical constant 22 22 USE eosbn2 ! Equation Of State 23 USE sbcmod ! ln_rnf 24 USE sbcrnf ! River runoff 23 USE sbcmod ! ln_rnf 24 USE sbcrnf ! River runoff 25 25 USE traqsr ! solar radiation penetration 26 26 USE trd_oce ! trends: ocean variables 27 USE trdtra ! trends manager: tracers 28 #if defined key_asminc 27 USE trdtra ! trends manager: tracers 28 #if defined key_asminc 29 29 USE asminc ! Assimilation increment 30 30 #endif … … 54 54 !!---------------------------------------------------------------------- 55 55 !! *** ROUTINE tra_sbc *** 56 !! 56 !! 57 57 !! ** Purpose : Compute the tracer surface boundary condition trend of 58 58 !! (flux through the interface, concentration/dilution effect) 59 59 !! and add it to the general trend of tracer equations. 60 60 !! 61 !! ** Method : The (air+ice)-sea flux has two components: 62 !! (1) Fext, external forcing (i.e. flux through the (air+ice)-sea interface); 63 !! (2) Fwe , tracer carried with the water that is exchanged with air+ice. 61 !! ** Method : The (air+ice)-sea flux has two components: 62 !! (1) Fext, external forcing (i.e. flux through the (air+ice)-sea interface); 63 !! (2) Fwe , tracer carried with the water that is exchanged with air+ice. 64 64 !! The input forcing fields (emp, rnf, sfx) contain Fext+Fwe, 65 65 !! they are simply added to the tracer trend (ts(Krhs)). … … 69 69 !! concentration/dilution effect associated with water exchanges. 70 70 !! 71 !! ** Action : - Update ts(Krhs) with the surface boundary condition trend 71 !! ** Action : - Update ts(Krhs) with the surface boundary condition trend 72 72 !! - send trends to trdtra module for further diagnostics(l_trdtra=T) 73 73 !!---------------------------------------------------------------------- … … 143 143 sbc_tsc(ji,jj,jp_sal) = r1_rho0 * sfx(ji,jj) ! salt flux due to freezing/melting 144 144 END_2D 145 IF( ln_linssh ) THEN !* linear free surface 145 IF( ln_linssh ) THEN !* linear free surface 146 146 DO_2D( isj, iej, isi, iei ) !==>> add concentration/dilution effect due to constant volume cell 147 147 sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm) … … 161 161 END_2D 162 162 END DO 163 ! 163 ! 164 164 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 165 165 IF( lrst_oce ) THEN !== write sbc_tsc in the ocean restart file ==! … … 173 173 !---------------------------------------- 174 174 ! 175 IF( ln_rnf ) THEN ! input of heat and salt due to river runoff 175 IF( ln_rnf ) THEN ! input of heat and salt due to river runoff 176 176 zfact = 0.5_wp 177 177 DO_2D( 0, 0, 0, 0 ) … … 182 182 & + ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep 183 183 IF( ln_rnf_sal ) pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) & 184 & + ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 184 & + ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 185 185 END DO 186 186 ENDIF … … 201 201 IF( ln_sshinc ) THEN ! input of heat and salt due to assimilation 202 202 ! 203 IF( ln_linssh ) THEN 203 IF( ln_linssh ) THEN 204 204 DO_2D( 0, 0, 0, 0 ) 205 205 ztim = ssh_iau(ji,jj) / e3t(ji,jj,1,Kmm) -
NEMO/trunk/src/OCE/TRA/trazdf.F90
r14010 r14072 17 17 USE phycst ! physical constant 18 18 USE zdf_oce ! ocean vertical physics variables 19 USE zdfmfc ! Mass FLux Convection 19 USE zdfmfc ! Mass FLux Convection 20 20 USE sbc_oce ! surface boundary condition: ocean 21 21 USE ldftra ! lateral diffusion: eddy diffusivity 22 USE ldfslp ! lateral diffusion: iso-neutral slope 22 USE ldfslp ! lateral diffusion: iso-neutral slope 23 23 USE trd_oce ! trends: ocean variables 24 24 USE trdtra ! trends: tracer trend manager … … 77 77 ! 78 78 ! !* compute lateral mixing trend and add it to the general trend 79 CALL tra_zdf_imp( kt, nit000, 'TRA', rDt, Kbb, Kmm, Krhs, pts, Kaa, jpts ) 79 CALL tra_zdf_imp( kt, nit000, 'TRA', rDt, Kbb, Kmm, Krhs, pts, Kaa, jpts ) 80 80 81 81 !!gm WHY here ! and I don't like that ! … … 113 113 END SUBROUTINE tra_zdf 114 114 115 116 SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt ) 115 116 SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt ) 117 117 !!---------------------------------------------------------------------- 118 118 !! *** ROUTINE tra_zdf_imp *** 119 119 !! 120 120 !! ** Purpose : Compute the after tracer through a implicit computation 121 !! of the vertical tracer diffusion (including the vertical component 122 !! of lateral mixing (only for 2nd order operator, for fourth order 123 !! it is already computed and add to the general trend in traldf) 121 !! of the vertical tracer diffusion (including the vertical component 122 !! of lateral mixing (only for 2nd order operator, for fourth order 123 !! it is already computed and add to the general trend in traldf) 124 124 !! 125 125 !! ** Method : The vertical diffusion of a tracer ,t , is given by: … … 169 169 zwt(:,:,1) = 0._wp 170 170 ! 171 IF( l_ldfslp ) THEN ! isoneutral diffusion: add the contribution 172 IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator 171 IF( l_ldfslp ) THEN ! isoneutral diffusion: add the contribution 172 IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator 173 173 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 174 zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk) 174 zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk) 175 175 END_3D 176 176 ELSE ! standard or triad iso-neutral operator … … 220 220 ! The solution will be in the 4d array pta. 221 221 ! The 3d array zwt is used as a work space array. 222 ! En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then 222 ! En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then 223 223 ! used as a work space array: its value is modified. 224 224 ! … … 230 230 END_3D 231 231 ! 232 ENDIF 233 ! 232 ENDIF 233 ! 234 234 ! Modification of rhs to add MF scheme 235 235 IF ( ln_zdfmfc ) THEN … … 239 239 DO_2D( 0, 0, 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 240 240 pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) & 241 & + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs) 241 & + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs) 242 242 END_2D 243 243 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) -
NEMO/trunk/src/OCE/TRA/zpshde.F90
r13982 r14072 7 7 !! NEMO 1.0 ! 2002-08 (G. Madec E. Durand) Optimization and Free form 8 8 !! - ! 2004-03 (C. Ethe) adapted for passive tracers 9 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA 9 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA 10 10 !! 3.6 ! 2014-11 (P. Mathiot) Add zps_hde_isf (needed to open a cavity) 11 11 !!====================================================================== 12 12 13 13 !!---------------------------------------------------------------------- 14 14 !! zps_hde : Horizontal DErivative of T, S and rd at the last … … 66 66 !!---------------------------------------------------------------------- 67 67 !! *** ROUTINE zps_hde *** 68 !! 68 !! 69 69 !! ** Purpose : Compute the horizontal derivative of T, S and rho 70 70 !! at u- and v-points with a linear interpolation for z-coordinate 71 71 !! with partial steps. 72 72 !! 73 !! ** Method : In z-coord with partial steps, scale factors on last 74 !! levels are different for each grid point, so that T, S and rd 73 !! ** Method : In z-coord with partial steps, scale factors on last 74 !! levels are different for each grid point, so that T, S and rd 75 75 !! points are not at the same depth as in z-coord. To have horizontal 76 !! gradients again, we interpolate T and S at the good depth : 77 !! Linear interpolation of T, S 76 !! gradients again, we interpolate T and S at the good depth : 77 !! Linear interpolation of T, S 78 78 !! Computation of di(tb) and dj(tb) by vertical interpolation: 79 79 !! di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~ 80 80 !! dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~ 81 81 !! This formulation computes the two cases: 82 !! CASE 1 CASE 2 82 !! CASE 1 CASE 2 83 83 !! k-1 ___ ___________ k-1 ___ ___________ 84 84 !! Ti T~ T~ Ti+1 … … 87 87 !! | |____ ____| | 88 88 !! ___ | | | ___ | | | 89 !! 89 !! 90 90 !! case 1-> e3w(i+1,:,:,Kmm) >= e3w(i,:,:,Kmm) ( and e3w(:,j+1,:,Kmm) >= e3w(:,j,:,Kmm) ) then 91 91 !! t~ = t(i+1,j ,k) + (e3w(i+1,j,k,Kmm) - e3w(i,j,k,Kmm)) * dk(Ti+1)/e3w(i+1,j,k,Kmm) … … 95 95 !! t~ = t(i,j,k) + (e3w(i,j,k,Kmm) - e3w(i+1,j,k,Kmm)) * dk(Ti)/e3w(i,j,k,Kmm) 96 96 !! ( t~ = t(i,j,k) + (e3w(i,j,k,Kmm) - e3w(i,j+1,k,Kmm)) * dk(Tj)/e3w(i,j,k,Kmm) ) 97 !! Idem for di(s) and dj(s) 97 !! Idem for di(s) and dj(s) 98 98 !! 99 99 !! For rho, we call eos which will compute rd~(t~,s~) at the right … … 175 175 ! 176 176 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp ) ! Lateral boundary cond. 177 ! 177 ! 178 178 IF( PRESENT( prd ) ) THEN !== horizontal derivative of density anomalies (rd) ==! (optional part) 179 179 pgru(:,:) = 0._wp … … 192 192 END_2D 193 193 ! 194 CALL eos( zti, zhi, zri ) ! interpolated density from zti, ztj 195 CALL eos( ztj, zhj, zrj ) ! at the partial step depth output in zri, zrj 194 CALL eos( zti, zhi, zri ) ! interpolated density from zti, ztj 195 CALL eos( ztj, zhj, zrj ) ! at the partial step depth output in zri, zrj 196 196 ! 197 197 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! Gradient of density at the last level … … 244 244 !!---------------------------------------------------------------------- 245 245 !! *** ROUTINE zps_hde_isf *** 246 !! 246 !! 247 247 !! ** Purpose : Compute the horizontal derivative of T, S and rho 248 248 !! at u- and v-points with a linear interpolation for z-coordinate 249 249 !! with partial steps for top (ice shelf) and bottom. 250 250 !! 251 !! ** Method : In z-coord with partial steps, scale factors on last 252 !! levels are different for each grid point, so that T, S and rd 251 !! ** Method : In z-coord with partial steps, scale factors on last 252 !! levels are different for each grid point, so that T, S and rd 253 253 !! points are not at the same depth as in z-coord. To have horizontal 254 254 !! gradients again, we interpolate T and S at the good depth : 255 255 !! For the bottom case: 256 !! Linear interpolation of T, S 256 !! Linear interpolation of T, S 257 257 !! Computation of di(tb) and dj(tb) by vertical interpolation: 258 258 !! di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~ 259 259 !! dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~ 260 260 !! This formulation computes the two cases: 261 !! CASE 1 CASE 2 261 !! CASE 1 CASE 2 262 262 !! k-1 ___ ___________ k-1 ___ ___________ 263 263 !! Ti T~ T~ Ti+1 … … 266 266 !! | |____ ____| | 267 267 !! ___ | | | ___ | | | 268 !! 268 !! 269 269 !! case 1-> e3w(i+1,j,k,Kmm) >= e3w(i,j,k,Kmm) ( and e3w(i,j+1,k,Kmm) >= e3w(i,j,k,Kmm) ) then 270 270 !! t~ = t(i+1,j ,k) + (e3w(i+1,j ,k,Kmm) - e3w(i,j,k,Kmm)) * dk(Ti+1)/e3w(i+1,j ,k,Kmm) … … 274 274 !! t~ = t(i,j,k) + (e3w(i,j,k,Kmm) - e3w(i+1,j ,k,Kmm)) * dk(Ti)/e3w(i,j,k,Kmm) 275 275 !! ( t~ = t(i,j,k) + (e3w(i,j,k,Kmm) - e3w(i ,j+1,k,Kmm)) * dk(Tj)/e3w(i,j,k,Kmm) ) 276 !! Idem for di(s) and dj(s) 276 !! Idem for di(s) and dj(s) 277 277 !! 278 278 !! For rho, we call eos which will compute rd~(t~,s~) at the right … … 364 364 ! horizontal derivative of density anomalies (rd) 365 365 IF( PRESENT( prd ) ) THEN ! depth of the partial step level 366 pgru(:,:)=0.0_wp ; pgrv(:,:)=0.0_wp ; 366 pgru(:,:)=0.0_wp ; pgrv(:,:)=0.0_wp ; 367 367 ! 368 368 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) … … 418 418 ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 419 419 ze3wu = gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm) 420 ze3wv = gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm) 420 ze3wv = gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm) 421 421 422 422 ! i- direction … … 463 463 ikv = mikv(ji,jj) 464 464 ze3wu = gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm) 465 ze3wv = gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm) 465 ze3wv = gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm) 466 466 ! 467 467 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = gdept(ji ,jj,iku,Kmm) ! i-direction: case 1 … … 475 475 END_2D 476 476 ! 477 CALL eos( zti, zhi, zri ) ! interpolated density from zti, ztj 478 CALL eos( ztj, zhj, zrj ) ! at the partial step depth output in zri, zrj 479 ! 480 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 481 iku = miku(ji,jj) 482 ikv = mikv(ji,jj) 477 CALL eos( zti, zhi, zri ) ! interpolated density from zti, ztj 478 CALL eos( ztj, zhj, zrj ) ! at the partial step depth output in zri, zrj 479 ! 480 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 481 iku = miku(ji,jj) 482 ikv = mikv(ji,jj) 483 483 ze3wu = gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm) 484 ze3wv = gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm) 484 ze3wv = gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm) 485 485 486 486 IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji ,jj ) - prd(ji,jj,iku) ) ! i: 1 … … 494 494 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'zpshde', pgrui, 'U', -1.0_wp , pgrvi, 'V', -1.0_wp ) ! Lateral boundary conditions 495 495 ! 496 END IF 496 END IF 497 497 ! 498 498 IF( ln_timing ) CALL timing_stop( 'zps_hde_isf')
Note: See TracChangeset
for help on using the changeset viewer.
