Changeset 12511 for NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA
- Timestamp:
- 2020-03-05T12:21:05+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/r12377_ticket2386
- Files:
-
- 13 edited
-
. (modified) (1 prop)
-
src/OCE/TRA/eosbn2.F90 (modified) (30 diffs)
-
src/OCE/TRA/traadv.F90 (modified) (2 diffs)
-
src/OCE/TRA/traadv_fct.F90 (modified) (1 diff)
-
src/OCE/TRA/traatf.F90 (modified) (9 diffs)
-
src/OCE/TRA/trabbc.F90 (modified) (4 diffs)
-
src/OCE/TRA/traldf_iso.F90 (modified) (3 diffs)
-
src/OCE/TRA/traldf_triad.F90 (modified) (3 diffs)
-
src/OCE/TRA/tramle.F90 (modified) (3 diffs)
-
src/OCE/TRA/tranpc.F90 (modified) (2 diffs)
-
src/OCE/TRA/traqsr.F90 (modified) (7 diffs)
-
src/OCE/TRA/trasbc.F90 (modified) (1 diff)
-
src/OCE/TRA/trazdf.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/r12377_ticket2386
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev _r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD ext/AGRIF5 ^/vendors/AGRIF/dev@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/eosbn2.F90
r12377 r12511 191 191 !! *** ROUTINE eos_insitu *** 192 192 !! 193 !! ** Purpose : Compute the in situ density (ratio rho/r au0) from193 !! ** Purpose : Compute the in situ density (ratio rho/rho0) from 194 194 !! potential temperature and salinity using an equation of state 195 195 !! selected in the nameos namelist 196 196 !! 197 !! ** Method : prd(t,s,z) = ( rho(t,s,z) - r au0 ) / rau0197 !! ** Method : prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0 198 198 !! with prd in situ density anomaly no units 199 199 !! t TEOS10: CT or EOS80: PT Celsius … … 201 201 !! z depth meters 202 202 !! rho in situ density kg/m^3 203 !! r au0 reference density kg/m^3203 !! rho0 reference density kg/m^3 204 204 !! 205 205 !! ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z). … … 210 210 !! 211 211 !! ln_seos : simplified equation of state 212 !! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / r au0212 !! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0 213 213 !! linear case function of T only: rn_alpha<>0, other coefficients = 0 214 214 !! linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0 … … 267 267 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 268 268 ! 269 prd(ji,jj,jk) = ( zn * r1_r au0 - 1._wp ) * ztm ! density anomaly (masked)269 prd(ji,jj,jk) = ( zn * r1_rho0 - 1._wp ) * ztm ! density anomaly (masked) 270 270 ! 271 271 END_3D … … 283 283 & - rn_nu * zt * zs 284 284 ! 285 prd(ji,jj,jk) = zn * r1_r au0 * ztm ! density anomaly (masked)285 prd(ji,jj,jk) = zn * r1_rho0 * ztm ! density anomaly (masked) 286 286 END_3D 287 287 ! … … 299 299 !! *** ROUTINE eos_insitu_pot *** 300 300 !! 301 !! ** Purpose : Compute the in situ density (ratio rho/r au0) and the301 !! ** Purpose : Compute the in situ density (ratio rho/rho0) and the 302 302 !! potential volumic mass (Kg/m3) from potential temperature and 303 303 !! salinity fields using an equation of state selected in the … … 379 379 prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp) ! potential density referenced at the surface 380 380 ! 381 prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_r au0 - 1._wp ) ! density anomaly (masked)381 prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_rho0 - 1._wp ) ! density anomaly (masked) 382 382 END DO 383 383 prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos … … 419 419 prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface 420 420 ! 421 prd(ji,jj,jk) = ( zn * r1_r au0 - 1._wp ) * ztm ! density anomaly (masked)421 prd(ji,jj,jk) = ( zn * r1_rho0 - 1._wp ) * ztm ! density anomaly (masked) 422 422 END_3D 423 423 ENDIF … … 434 434 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs & 435 435 & - rn_nu * zt * zs 436 prhop(ji,jj,jk) = ( r au0 + zn ) * ztm436 prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 437 437 ! ! density anomaly (masked) 438 438 zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 439 prd(ji,jj,jk) = zn * r1_r au0 * ztm439 prd(ji,jj,jk) = zn * r1_rho0 * ztm 440 440 ! 441 441 END_3D … … 454 454 !! *** ROUTINE eos_insitu_2d *** 455 455 !! 456 !! ** Purpose : Compute the in situ density (ratio rho/r au0) from456 !! ** Purpose : Compute the in situ density (ratio rho/rho0) from 457 457 !! potential temperature and salinity using an equation of state 458 458 !! selected in the nameos namelist. * 2D field case … … 508 508 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 509 509 ! 510 prd(ji,jj) = zn * r1_r au0 - 1._wp ! unmasked in situ density anomaly510 prd(ji,jj) = zn * r1_rho0 - 1._wp ! unmasked in situ density anomaly 511 511 ! 512 512 END_2D … … 524 524 & - rn_nu * zt * zs 525 525 ! 526 prd(ji,jj) = zn * r1_r au0 ! unmasked in situ density anomaly526 prd(ji,jj) = zn * r1_rho0 ! unmasked in situ density anomaly 527 527 ! 528 528 END_2D … … 588 588 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 589 589 ! 590 pab(ji,jj,jk,jp_tem) = zn * r1_r au0 * ztm590 pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm 591 591 ! 592 592 ! beta … … 609 609 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 610 610 ! 611 pab(ji,jj,jk,jp_sal) = zn / zs * r1_r au0 * ztm611 pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm 612 612 ! 613 613 END_3D … … 622 622 ! 623 623 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 624 pab(ji,jj,jk,jp_tem) = zn * r1_r au0 * ztm ! alpha624 pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm ! alpha 625 625 ! 626 626 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 627 pab(ji,jj,jk,jp_sal) = zn * r1_r au0 * ztm ! beta627 pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm ! beta 628 628 ! 629 629 END_3D … … 694 694 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 695 695 ! 696 pab(ji,jj,jp_tem) = zn * r1_r au0696 pab(ji,jj,jp_tem) = zn * r1_rho0 697 697 ! 698 698 ! beta … … 715 715 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 716 716 ! 717 pab(ji,jj,jp_sal) = zn / zs * r1_r au0717 pab(ji,jj,jp_sal) = zn / zs * r1_rho0 718 718 ! 719 719 ! … … 729 729 ! 730 730 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 731 pab(ji,jj,jp_tem) = zn * r1_r au0 ! alpha731 pab(ji,jj,jp_tem) = zn * r1_rho0 ! alpha 732 732 ! 733 733 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 734 pab(ji,jj,jp_sal) = zn * r1_r au0 ! beta734 pab(ji,jj,jp_sal) = zn * r1_rho0 ! beta 735 735 ! 736 736 END_2D … … 799 799 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 800 800 ! 801 pab(jp_tem) = zn * r1_r au0801 pab(jp_tem) = zn * r1_rho0 802 802 ! 803 803 ! beta … … 820 820 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 821 821 ! 822 pab(jp_sal) = zn / zs * r1_r au0822 pab(jp_sal) = zn / zs * r1_rho0 823 823 ! 824 824 ! … … 831 831 ! 832 832 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 833 pab(jp_tem) = zn * r1_r au0 ! alpha833 pab(jp_tem) = zn * r1_rho0 ! alpha 834 834 ! 835 835 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 836 pab(jp_sal) = zn * r1_r au0 ! beta836 pab(jp_sal) = zn * r1_rho0 ! beta 837 837 ! 838 838 CASE DEFAULT … … 1052 1052 !! ** Method : PE is defined analytically as the vertical 1053 1053 !! primitive of EOS times -g integrated between 0 and z>0. 1054 !! pen is the nonlinear bsq-PE anomaly: pen = ( PE - r au0 gz ) / rau0 gz - rd1054 !! pen is the nonlinear bsq-PE anomaly: pen = ( PE - rho0 gz ) / rho0 gz - rd 1055 1055 !! = 1/z * /int_0^z rd dz - rd 1056 1056 !! where rd is the density anomaly (see eos_rhd function) 1057 1057 !! ab_pe are partial derivatives of PE anomaly with respect to T and S: 1058 !! ab_pe(1) = - 1/(r au0 gz) * dPE/dT + drd/dT = - d(pen)/dT1059 !! ab_pe(2) = 1/(r au0 gz) * dPE/dS + drd/dS = d(pen)/dS1058 !! ab_pe(1) = - 1/(rho0 gz) * dPE/dT + drd/dT = - d(pen)/dT 1059 !! ab_pe(2) = 1/(rho0 gz) * dPE/dS + drd/dS = d(pen)/dS 1060 1060 !! 1061 1061 !! ** Action : - pen : PE anomaly given at T-points … … 1103 1103 zn = ( zn2 * zh + zn1 ) * zh + zn0 1104 1104 ! 1105 ppen(ji,jj,jk) = zn * zh * r1_r au0 * ztm1105 ppen(ji,jj,jk) = zn * zh * r1_rho0 * ztm 1106 1106 ! 1107 1107 ! alphaPE non-linear anomaly … … 1118 1118 zn = ( zn2 * zh + zn1 ) * zh + zn0 1119 1119 ! 1120 pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_r au0 * ztm1120 pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm 1121 1121 ! 1122 1122 ! betaPE non-linear anomaly … … 1133 1133 zn = ( zn2 * zh + zn1 ) * zh + zn0 1134 1134 ! 1135 pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_r au0 * ztm1135 pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm 1136 1136 ! 1137 1137 END_3D … … 1144 1144 zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point 1145 1145 ztm = tmask(ji,jj,jk) ! tmask 1146 zn = 0.5_wp * zh * r1_r au0 * ztm1146 zn = 0.5_wp * zh * r1_rho0 * ztm 1147 1147 ! ! Potential Energy 1148 1148 ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn … … 1186 1186 IF(lwm) WRITE( numond, nameos ) 1187 1187 ! 1188 r au0 = 1026._wp !: volumic mass of reference [kg/m3]1188 rho0 = 1026._wp !: volumic mass of reference [kg/m3] 1189 1189 rcp = 3991.86795711963_wp !: heat capacity [J/K] 1190 1190 ! … … 1598 1598 WRITE(numout,*) ' ==>>> use of simplified eos: ' 1599 1599 WRITE(numout,*) ' rhd(dT=T-10,dS=S-35,Z) = [-a0*(1+lambda1/2*dT+mu1*Z)*dT ' 1600 WRITE(numout,*) ' + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / r au0'1600 WRITE(numout,*) ' + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / rho0' 1601 1601 WRITE(numout,*) ' with the following coefficients :' 1602 1602 WRITE(numout,*) ' thermal exp. coef. rn_a0 = ', rn_a0 … … 1617 1617 END SELECT 1618 1618 ! 1619 r au0_rcp = rau0 * rcp1620 r1_r au0 = 1._wp / rau01619 rho0_rcp = rho0 * rcp 1620 r1_rho0 = 1._wp / rho0 1621 1621 r1_rcp = 1._wp / rcp 1622 r1_r au0_rcp = 1._wp / rau0_rcp1622 r1_rho0_rcp = 1._wp / rho0_rcp 1623 1623 ! 1624 1624 IF(lwp) THEN … … 1635 1635 IF(lwp) WRITE(numout,*) 1636 1636 IF(lwp) WRITE(numout,*) ' Associated physical constant' 1637 IF(lwp) WRITE(numout,*) ' volumic mass of reference r au0 = ', rau0 , ' kg/m^3'1638 IF(lwp) WRITE(numout,*) ' 1. / r au0 r1_rau0 = ', r1_rau0, ' m^3/kg'1637 IF(lwp) WRITE(numout,*) ' volumic mass of reference rho0 = ', rho0 , ' kg/m^3' 1638 IF(lwp) WRITE(numout,*) ' 1. / rho0 r1_rho0 = ', r1_rho0, ' m^3/kg' 1639 1639 IF(lwp) WRITE(numout,*) ' ocean specific heat rcp = ', rcp , ' J/Kelvin' 1640 IF(lwp) WRITE(numout,*) ' r au0 * rcp rau0_rcp = ', rau0_rcp1641 IF(lwp) WRITE(numout,*) ' 1. / ( r au0 * rcp ) r1_rau0_rcp = ', r1_rau0_rcp1640 IF(lwp) WRITE(numout,*) ' rho0 * rcp rho0_rcp = ', rho0_rcp 1641 IF(lwp) WRITE(numout,*) ' 1. / ( rho0 * rcp ) r1_rho0_rcp = ', r1_rho0_rcp 1642 1642 ! 1643 1643 END SUBROUTINE eos_init -
NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traadv.F90
r12377 r12511 92 92 IF( ln_timing ) CALL timing_start('tra_adv') 93 93 ! 94 ! ! set time step95 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! at nit000 (Euler)96 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2._wp * rdt ! at nit000 or nit000+1 (Leapfrog)97 ENDIF98 !99 94 ! !== effective transport ==! 100 95 zuu(:,:,jpk) = 0._wp … … 149 144 CALL tra_adv_cen ( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v ) 150 145 CASE ( np_FCT ) ! FCT scheme : 2nd / 4th order 151 CALL tra_adv_fct ( kt, nit000, 'TRA', r 2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v )146 CALL tra_adv_fct ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v ) 152 147 CASE ( np_MUS ) ! MUSCL 153 CALL tra_adv_mus ( kt, nit000, 'TRA', r 2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )148 CALL tra_adv_mus ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 154 149 CASE ( np_UBS ) ! UBS 155 CALL tra_adv_ubs ( kt, nit000, 'TRA', r 2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v )150 CALL tra_adv_ubs ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v ) 156 151 CASE ( np_QCK ) ! QUICKEST 157 CALL tra_adv_qck ( kt, nit000, 'TRA', r 2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs )152 CALL tra_adv_qck ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs ) 158 153 ! 159 154 END SELECT -
NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traadv_fct.F90
r12377 r12511 20 20 USE diaptr ! poleward transport diagnostics 21 21 USE diaar5 ! AR5 diagnostics 22 USE phycst , ONLY : r au0_rcp22 USE phycst , ONLY : rho0_rcp 23 23 USE zdf_oce , ONLY : ln_zad_Aimp 24 24 ! -
NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traatf.F90
r12377 r12511 113 113 IF( ln_bdy ) CALL bdy_tra( kt, Kbb, pts, Kaa ) ! BDY open boundaries 114 114 115 ! set time step size (Euler/Leapfrog)116 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! at nit000 (Euler)117 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2._wp* rdt ! at nit000 or nit000+1 (Leapfrog)118 ENDIF119 120 115 ! trends computation initialisation 121 116 IF( l_trdtra ) THEN … … 128 123 ENDIF 129 124 ! total trend for the non-time-filtered variables. 130 zfact = 1.0 / r dt125 zfact = 1.0 / rn_Dt 131 126 ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from pts(Kmm) terms 132 127 DO jk = 1, jpkm1 … … 144 139 ENDIF 145 140 146 IF( neuler == 0 .AND. kt == nit000) THEN ! Euler time-stepping141 IF( l_1st_euler ) THEN ! Euler time-stepping 147 142 ! 148 143 IF (l_trdtra .AND. .NOT. ln_linssh ) THEN ! Zero Asselin filter contribution must be explicitly written out since for vvl … … 156 151 ELSE ! Leap-Frog + Asselin filter time stepping 157 152 ! 158 IF( ln_linssh ) THEN ; CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000, 'TRA', pts, jpts ) ! linear free surface159 ELSE ; CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nit000, r dt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts ) ! non-linear free surface153 IF( ln_linssh ) THEN ; CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000, 'TRA', pts, jpts ) ! linear free surface 154 ELSE ; CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nit000, rn_Dt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts ) ! non-linear free surface 160 155 ENDIF 161 156 ! … … 167 162 ! 168 163 IF( l_trdtra .AND. ln_linssh ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt 169 zfact = 1._wp / r 2dt164 zfact = 1._wp / rDt 170 165 DO jk = 1, jpkm1 171 166 ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kmm) - ztrdt(:,:,jk) ) * zfact … … 219 214 ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb) ! time laplacian on tracers 220 215 ! 221 pt(ji,jj,jk,jn,Kmm) = ztn + atfp * ztd ! pt <-- filtered pt216 pt(ji,jj,jk,jn,Kmm) = ztn + rn_atfp * ztd ! pt <-- filtered pt 222 217 END_3D 223 218 ! … … 234 229 !! 235 230 !! ** Method : - Apply a thickness weighted Asselin time filter on now fields. 236 !! pt(Kmm) = ( e3t(Kmm)*pt(Kmm) + atfp*[ e3t(Kbb)*pt(Kbb) - 2 e3t(Kmm)*pt(Kmm) + e3t_a*pt(Kaa) ] )237 !! /( e3t(Kmm) + atfp*[ e3t(Kbb) - 2 e3t(Kmm) + e3t(Kaa) ] )231 !! pt(Kmm) = ( e3t(Kmm)*pt(Kmm) + rn_atfp*[ e3t(Kbb)*pt(Kbb) - 2 e3t(Kmm)*pt(Kmm) + e3t_a*pt(Kaa) ] ) 232 !! /( e3t(Kmm) + rn_atfp*[ e3t(Kbb) - 2 e3t(Kmm) + e3t(Kaa) ] ) 238 233 !! 239 234 !! ** Action : - pt(Kmm) ready for the next time step … … 277 272 ENDIF 278 273 zfact = 1._wp / p2dt 279 zfact1 = atfp * p2dt280 zfact2 = zfact1 * r1_r au0274 zfact1 = rn_atfp * p2dt 275 zfact2 = zfact1 * r1_rho0 281 276 DO jn = 1, kjpt 282 277 DO_3D_00_00( 1, jpkm1 ) … … 292 287 ztc_d = ztc_a - 2. * ztc_n + ztc_b 293 288 ! 294 ze3t_f = ze3t_n + atfp * ze3t_d295 ztc_f = ztc_n + atfp * ztc_d289 ze3t_f = ze3t_n + rn_atfp * ze3t_d 290 ztc_f = ztc_n + rn_atfp * ztc_d 296 291 ! 297 292 ! Add asselin correction on scale factors: -
NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/trabbc.F90
r12377 r12511 66 66 !! ocean bottom can be computed once and is added to the temperature 67 67 !! trend juste above the bottom at each time step: 68 !! ta = ta + Qsf / (r au0 rcp e3T) for k= mbkt68 !! ta = ta + Qsf / (rho0 rcp e3T) for k= mbkt 69 69 !! Where Qsf is the geothermal heat flux. 70 70 !! … … 102 102 ENDIF 103 103 ! 104 CALL iom_put ( "hfgeou" , r au0_rcp * qgh_trd0(:,:) )104 CALL iom_put ( "hfgeou" , rho0_rcp * qgh_trd0(:,:) ) 105 105 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbc - Ta: ', mask1=tmask, clinfo3='tra-ta' ) 106 106 ! … … 162 162 CASE ( 1 ) !* constant flux 163 163 IF(lwp) WRITE(numout,*) ' ==>>> constant heat flux = ', rn_geoflx_cst 164 qgh_trd0(:,:) = r1_r au0_rcp * rn_geoflx_cst164 qgh_trd0(:,:) = r1_rho0_rcp * rn_geoflx_cst 165 165 ! 166 166 CASE ( 2 ) !* variable geothermal heat flux : read the geothermal fluxes in mW/m2 … … 179 179 180 180 CALL fld_read( nit000, 1, sf_qgh ) ! Read qgh data 181 qgh_trd0(:,:) = r1_r au0_rcp * sf_qgh(1)%fnow(:,:,1) * 1.e-3 ! conversion in W/m2181 qgh_trd0(:,:) = r1_rho0_rcp * sf_qgh(1)%fnow(:,:,1) * 1.e-3 ! conversion in W/m2 182 182 ! 183 183 CASE DEFAULT -
NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traldf_iso.F90
r12377 r12511 109 109 REAL(wp) :: zmsku, zahu_w, zabe1, zcof1, zcoef3 ! local scalars 110 110 REAL(wp) :: zmskv, zahv_w, zabe2, zcof2, zcoef4 ! - - 111 REAL(wp) :: zcoef0, ze3w_2, zsign , z2dt, z1_2dt! - -111 REAL(wp) :: zcoef0, ze3w_2, zsign ! - - 112 112 REAL(wp), DIMENSION(jpi,jpj) :: zdkt, zdk1t, z2d 113 113 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdit, zdjt, zftu, zftv, ztfw … … 129 129 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 130 130 ! 131 ! ! set time step size (Euler/Leapfrog)132 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2dt = rdt ! at nit000 (Euler)133 ELSE ; z2dt = 2.* rdt ! (Leapfrog)134 ENDIF135 z1_2dt = 1._wp / z2dt136 131 ! 137 132 IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0) … … 178 173 DO_3D_10_10( 2, jpkm1 ) 179 174 ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 180 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 )181 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt175 zcoef0 = rDt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 176 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_Dt 182 177 END_3D 183 178 ENDIF -
NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traldf_triad.F90
r12377 r12511 86 86 INTEGER :: ip,jp,kp ! dummy loop indices 87 87 INTEGER :: ierr ! local integer 88 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 89 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 90 REAL(wp) :: zcoef0, ze3w_2, zsign , z2dt, z1_2dt! - -88 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars 89 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 90 REAL(wp) :: zcoef0, ze3w_2, zsign ! - - 91 91 ! 92 92 REAL(wp) :: zslope_skew, zslope_iso, zslope2, zbu, zbv … … 111 111 l_hst = .FALSE. 112 112 l_ptr = .FALSE. 113 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) ) l_ptr = .TRUE. 114 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 115 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 116 ! 117 ! ! set time step size (Euler/Leapfrog) 118 IF( neuler == 0 .AND. kt == kit000 ) THEN ; z2dt = rdt ! at nit000 (Euler) 119 ELSE ; z2dt = 2.* rdt ! (Leapfrog) 113 IF( cdtype == 'TRA' ) THEN 114 IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf') ) l_ptr = .TRUE. 115 IF( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 116 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) l_hst = .TRUE. 120 117 ENDIF 121 z1_2dt = 1._wp / z2dt122 118 ! 123 119 IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0) … … 189 185 DO_3D_10_10( 2, jpkm1 ) 190 186 ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 191 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 )192 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt187 zcoef0 = rDt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 188 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_Dt 193 189 END_3D 194 190 ENDIF -
NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/tramle.F90
r12377 r12511 41 41 42 42 REAL(wp) :: r5_21 = 5.e0 / 21.e0 ! factor used in mle streamfunction computation 43 REAL(wp) :: rb_c ! ML buoyancy criteria = g rho_c /r au0 where rho_c is defined in zdfmld43 REAL(wp) :: rb_c ! ML buoyancy criteria = g rho_c /rho0 where rho_c is defined in zdfmld 44 44 REAL(wp) :: rc_f ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_mle=1 case 45 45 … … 112 112 zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 113 113 zmld(ji,jj) = zmld(ji,jj) + zc 114 zbm (ji,jj) = zbm (ji,jj) + zc * (r au0 - rhop(ji,jj,jk) ) * r1_rau0114 zbm (ji,jj) = zbm (ji,jj) + zc * (rho0 - rhop(ji,jj,jk) ) * r1_rho0 115 115 zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 116 116 END_3D … … 273 273 IF( ln_mle ) THEN ! MLE initialisation 274 274 ! 275 rb_c = grav * rn_rho_c_mle /r au0 ! Mixed Layer buoyancy criteria275 rb_c = grav * rn_rho_c_mle /rho0 ! Mixed Layer buoyancy criteria 276 276 IF(lwp) WRITE(numout,*) 277 277 IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 ' -
NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/tranpc.F90
r12377 r12511 67 67 LOGICAL :: l_bottom_reached, l_column_treated 68 68 REAL(wp) :: zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z 69 REAL(wp) :: zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r 2dt69 REAL(wp) :: zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_rDt 70 70 REAL(wp), PARAMETER :: zn2_zero = 1.e-14_wp ! acceptance criteria for neutrality (N2==0) 71 71 REAL(wp), DIMENSION( jpk ) :: zvn2 ! vertical profile of N2 at 1 given point... … … 301 301 ! 302 302 IF( l_trdtra ) THEN ! send the Non penetrative mixing trends for diagnostic 303 z1_r 2dt = 1._wp / (2._wp * rdt)304 ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * z1_r 2dt305 ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * z1_r 2dt303 z1_rDt = 1._wp / (2._wp * rn_Dt) 304 ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * z1_rDt 305 ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * z1_rDt 306 306 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_npc, ztrdt ) 307 307 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_npc, ztrds ) -
NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traqsr.F90
r12377 r12511 87 87 !! I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 88 88 !! The temperature trend associated with the solar radiation penetration 89 !! is given by : zta = 1/e3t dk[ I ] / (r au0*Cp)89 !! is given by : zta = 1/e3t dk[ I ] / (rho0*Cp) 90 90 !! At the bottom, boudary condition for the radiation is no flux : 91 91 !! all heat which has not been absorbed in the above levels is put … … 135 135 ! !-----------------------------------! 136 136 IF( kt == nit000 ) THEN !== 1st time step ==! 137 !!gm case neuler not taken into account.... 138 IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN ! read in restart 137 IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 .AND. .NOT.l_1st_euler ) THEN ! read in restart 139 138 IF(lwp) WRITE(numout,*) ' nit000-1 qsr tracer content forcing field read in the restart file' 140 139 z1_2 = 0.5_wp … … 156 155 ! 157 156 DO jk = 1, nksr 158 qsr_hc(:,:,jk) = r1_r au0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )157 qsr_hc(:,:,jk) = r1_rho0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) ) 159 158 END DO 160 159 ! … … 228 227 ! 229 228 DO_3D_00_00( 1, nksr ) 230 qsr_hc(ji,jj,jk) = r1_r au0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) )229 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) ) 231 230 END_3D 232 231 ! … … 235 234 CASE( np_2BD ) !== 2-bands fluxes ==! 236 235 ! 237 zz0 = rn_abs * r1_r au0_rcp ! surface equi-partition in 2-bands238 zz1 = ( 1. - rn_abs ) * r1_r au0_rcp236 zz0 = rn_abs * r1_rho0_rcp ! surface equi-partition in 2-bands 237 zz1 = ( 1. - rn_abs ) * r1_rho0_rcp 239 238 DO_3D_00_00( 1, nksr ) 240 239 zc0 = zz0 * EXP( -gdepw(ji,jj,jk ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk ,Kmm)*xsi1r ) … … 253 252 ! sea-ice: store the 1st ocean level attenuation coefficient 254 253 DO_2D_00_00 255 IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_r au0_rcp * qsr(ji,jj) )254 IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) ) 256 255 ELSE ; fraqsr_1lev(ji,jj) = 1._wp 257 256 ENDIF … … 263 262 zetot(:,:,nksr+1:jpk) = 0._wp ! below ~400m set to zero 264 263 DO jk = nksr, 1, -1 265 zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * r au0_rcp264 zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp 266 265 END DO 267 266 CALL iom_put( 'qsr3d', zetot ) ! 3D distribution of shortwave Radiation -
NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/trasbc.F90
r12377 r12511 124 124 ! !== Now sbc tracer content fields ==! 125 125 DO_2D_01_00 126 sbc_tsc(ji,jj,jp_tem) = r1_r au0_rcp * qns(ji,jj) ! non solar heat flux127 sbc_tsc(ji,jj,jp_sal) = r1_r au0 * sfx(ji,jj) ! salt flux due to freezing/melting126 sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj) ! non solar heat flux 127 sbc_tsc(ji,jj,jp_sal) = r1_rho0 * sfx(ji,jj) ! salt flux due to freezing/melting 128 128 END_2D 129 129 IF( ln_linssh ) THEN !* linear free surface 130 130 DO_2D_01_00 131 sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_r au0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm)132 sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_r au0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm)131 sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm) 132 sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm) 133 133 END_2D 134 134 IF( iom_use('emp_x_sst') ) CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) ) -
NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/trazdf.F90
r12377 r12511 66 66 ENDIF 67 67 ! 68 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! at nit000, = rdt (restarting with Euler time stepping)69 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2. * rdt ! otherwise, = 2 rdt (leapfrog)70 ENDIF71 !72 68 IF( l_trdtra ) THEN !* Save ta and sa trends 73 69 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) … … 77 73 ! 78 74 ! !* compute lateral mixing trend and add it to the general trend 79 CALL tra_zdf_imp( kt, nit000, 'TRA', r 2dt, Kbb, Kmm, Krhs, pts, Kaa, jpts )75 CALL tra_zdf_imp( kt, nit000, 'TRA', rDt, Kbb, Kmm, Krhs, pts, Kaa, jpts ) 80 76 81 77 !!gm WHY here ! and I don't like that ! … … 89 85 DO jk = 1, jpkm1 90 86 ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) & 91 & / (e3t(:,:,jk,Kmm)*r 2dt) ) - ztrdt(:,:,jk)87 & / (e3t(:,:,jk,Kmm)*rDt) ) - ztrdt(:,:,jk) 92 88 ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) & 93 & / (e3t(:,:,jk,Kmm)*r 2dt) ) - ztrds(:,:,jk)89 & / (e3t(:,:,jk,Kmm)*rDt) ) - ztrds(:,:,jk) 94 90 END DO 95 91 !!gm this should be moved in trdtra.F90 and done on all trends
Note: See TracChangeset
for help on using the changeset viewer.
