Changeset 9939 for NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA
- Timestamp:
- 2018-07-13T09:28:50+02:00 (6 years ago)
- Location:
- NEMO/branches/2018/dev_r9838_ENHANCE04_RK3
- Files:
-
- 15 edited
- 1 copied
-
. (copied) (copied from NEMO/trunk)
-
src/OCE/TRA/eosbn2.F90 (modified) (30 diffs)
-
src/OCE/TRA/traadv.F90 (modified) (4 diffs)
-
src/OCE/TRA/traadv_fct.F90 (modified) (3 diffs)
-
src/OCE/TRA/trabbc.F90 (modified) (6 diffs)
-
src/OCE/TRA/trabbl.F90 (modified) (3 diffs)
-
src/OCE/TRA/tradmp.F90 (modified) (3 diffs)
-
src/OCE/TRA/traldf.F90 (modified) (3 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) (4 diffs)
-
src/OCE/TRA/tranxt.F90 (modified) (9 diffs)
-
src/OCE/TRA/traqsr.F90 (modified) (10 diffs)
-
src/OCE/TRA/trasbc.F90 (modified) (7 diffs)
-
src/OCE/TRA/trazdf.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/eosbn2.F90
r9757 r9939 190 190 !! *** ROUTINE eos_insitu *** 191 191 !! 192 !! ** Purpose : Compute the in situ density (ratio rho/r au0) from192 !! ** Purpose : Compute the in situ density (ratio rho/rho0) from 193 193 !! potential temperature and salinity using an equation of state 194 194 !! selected in the nameos namelist 195 195 !! 196 !! ** Method : prd(t,s,z) = ( rho(t,s,z) - r au0 ) / rau0196 !! ** Method : prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0 197 197 !! with prd in situ density anomaly no units 198 198 !! t TEOS10: CT or EOS80: PT Celsius … … 200 200 !! z depth meters 201 201 !! rho in situ density kg/m^3 202 !! r au0 reference density kg/m^3202 !! rho0 reference density kg/m^3 203 203 !! 204 204 !! ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z). … … 209 209 !! 210 210 !! ln_seos : simplified equation of state 211 !! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / r au0211 !! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0 212 212 !! linear case function of T only: rn_alpha<>0, other coefficients = 0 213 213 !! linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0 … … 268 268 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 269 269 ! 270 prd(ji,jj,jk) = ( zn * r1_r au0 - 1._wp ) * ztm ! density anomaly (masked)270 prd(ji,jj,jk) = ( zn * r1_rho0 - 1._wp ) * ztm ! density anomaly (masked) 271 271 ! 272 272 END DO … … 288 288 & - rn_nu * zt * zs 289 289 ! 290 prd(ji,jj,jk) = zn * r1_r au0 * ztm ! density anomaly (masked)290 prd(ji,jj,jk) = zn * r1_rho0 * ztm ! density anomaly (masked) 291 291 END DO 292 292 END DO … … 306 306 !! *** ROUTINE eos_insitu_pot *** 307 307 !! 308 !! ** Purpose : Compute the in situ density (ratio rho/r au0) and the308 !! ** Purpose : Compute the in situ density (ratio rho/rho0) and the 309 309 !! potential volumic mass (Kg/m3) from potential temperature and 310 310 !! salinity fields using an equation of state selected in the … … 388 388 prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp) ! potential density referenced at the surface 389 389 ! 390 prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_r au0 - 1._wp ) ! density anomaly (masked)390 prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_rho0 - 1._wp ) ! density anomaly (masked) 391 391 END DO 392 392 prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos … … 432 432 prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface 433 433 ! 434 prd(ji,jj,jk) = ( zn * r1_r au0 - 1._wp ) * ztm ! density anomaly (masked)434 prd(ji,jj,jk) = ( zn * r1_rho0 - 1._wp ) * ztm ! density anomaly (masked) 435 435 END DO 436 436 END DO … … 451 451 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs & 452 452 & - rn_nu * zt * zs 453 prhop(ji,jj,jk) = ( r au0 + zn ) * ztm453 prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 454 454 ! ! density anomaly (masked) 455 455 zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 456 prd(ji,jj,jk) = zn * r1_r au0 * ztm456 prd(ji,jj,jk) = zn * r1_rho0 * ztm 457 457 ! 458 458 END DO … … 473 473 !! *** ROUTINE eos_insitu_2d *** 474 474 !! 475 !! ** Purpose : Compute the in situ density (ratio rho/r au0) from475 !! ** Purpose : Compute the in situ density (ratio rho/rho0) from 476 476 !! potential temperature and salinity using an equation of state 477 477 !! selected in the nameos namelist. * 2D field case … … 528 528 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 529 529 ! 530 prd(ji,jj) = zn * r1_r au0 - 1._wp ! unmasked in situ density anomaly530 prd(ji,jj) = zn * r1_rho0 - 1._wp ! unmasked in situ density anomaly 531 531 ! 532 532 END DO … … 548 548 & - rn_nu * zt * zs 549 549 ! 550 prd(ji,jj) = zn * r1_r au0 ! unmasked in situ density anomaly550 prd(ji,jj) = zn * r1_rho0 ! unmasked in situ density anomaly 551 551 ! 552 552 END DO … … 616 616 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 617 617 ! 618 pab(ji,jj,jk,jp_tem) = zn * r1_r au0 * ztm618 pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm 619 619 ! 620 620 ! beta … … 637 637 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 638 638 ! 639 pab(ji,jj,jk,jp_sal) = zn / zs * r1_r au0 * ztm639 pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm 640 640 ! 641 641 END DO … … 654 654 ! 655 655 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 656 pab(ji,jj,jk,jp_tem) = zn * r1_r au0 * ztm ! alpha656 pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm ! alpha 657 657 ! 658 658 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 659 pab(ji,jj,jk,jp_sal) = zn * r1_r au0 * ztm ! beta659 pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm ! beta 660 660 ! 661 661 END DO … … 729 729 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 730 730 ! 731 pab(ji,jj,jp_tem) = zn * r1_r au0731 pab(ji,jj,jp_tem) = zn * r1_rho0 732 732 ! 733 733 ! beta … … 750 750 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 751 751 ! 752 pab(ji,jj,jp_sal) = zn / zs * r1_r au0752 pab(ji,jj,jp_sal) = zn / zs * r1_rho0 753 753 ! 754 754 ! … … 768 768 ! 769 769 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 770 pab(ji,jj,jp_tem) = zn * r1_r au0 ! alpha770 pab(ji,jj,jp_tem) = zn * r1_rho0 ! alpha 771 771 ! 772 772 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 773 pab(ji,jj,jp_sal) = zn * r1_r au0 ! beta773 pab(ji,jj,jp_sal) = zn * r1_rho0 ! beta 774 774 ! 775 775 END DO … … 841 841 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 842 842 ! 843 pab(jp_tem) = zn * r1_r au0843 pab(jp_tem) = zn * r1_rho0 844 844 ! 845 845 ! beta … … 862 862 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 863 863 ! 864 pab(jp_sal) = zn / zs * r1_r au0864 pab(jp_sal) = zn / zs * r1_rho0 865 865 ! 866 866 ! … … 873 873 ! 874 874 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 875 pab(jp_tem) = zn * r1_r au0 ! alpha875 pab(jp_tem) = zn * r1_rho0 ! alpha 876 876 ! 877 877 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 878 pab(jp_sal) = zn * r1_r au0 ! beta878 pab(jp_sal) = zn * r1_rho0 ! beta 879 879 ! 880 880 CASE DEFAULT … … 1104 1104 !! ** Method : PE is defined analytically as the vertical 1105 1105 !! primitive of EOS times -g integrated between 0 and z>0. 1106 !! pen is the nonlinear bsq-PE anomaly: pen = ( PE - r au0 gz ) / rau0 gz - rd1106 !! pen is the nonlinear bsq-PE anomaly: pen = ( PE - rho0 gz ) / rho0 gz - rd 1107 1107 !! = 1/z * /int_0^z rd dz - rd 1108 1108 !! where rd is the density anomaly (see eos_rhd function) 1109 1109 !! ab_pe are partial derivatives of PE anomaly with respect to T and S: 1110 !! ab_pe(1) = - 1/(r au0 gz) * dPE/dT + drd/dT = - d(pen)/dT1111 !! ab_pe(2) = 1/(r au0 gz) * dPE/dS + drd/dS = d(pen)/dS1110 !! ab_pe(1) = - 1/(rho0 gz) * dPE/dT + drd/dT = - d(pen)/dT 1111 !! ab_pe(2) = 1/(rho0 gz) * dPE/dS + drd/dS = d(pen)/dS 1112 1112 !! 1113 1113 !! ** Action : - pen : PE anomaly given at T-points … … 1156 1156 zn = ( zn2 * zh + zn1 ) * zh + zn0 1157 1157 ! 1158 ppen(ji,jj,jk) = zn * zh * r1_r au0 * ztm1158 ppen(ji,jj,jk) = zn * zh * r1_rho0 * ztm 1159 1159 ! 1160 1160 ! alphaPE non-linear anomaly … … 1171 1171 zn = ( zn2 * zh + zn1 ) * zh + zn0 1172 1172 ! 1173 pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_r au0 * ztm1173 pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm 1174 1174 ! 1175 1175 ! betaPE non-linear anomaly … … 1186 1186 zn = ( zn2 * zh + zn1 ) * zh + zn0 1187 1187 ! 1188 pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_r au0 * ztm1188 pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm 1189 1189 ! 1190 1190 END DO … … 1201 1201 zh = gdept_n(ji,jj,jk) ! depth in meters at t-point 1202 1202 ztm = tmask(ji,jj,jk) ! tmask 1203 zn = 0.5_wp * zh * r1_r au0 * ztm1203 zn = 0.5_wp * zh * r1_rho0 * ztm 1204 1204 ! ! Potential Energy 1205 1205 ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn … … 1248 1248 IF(lwm) WRITE( numond, nameos ) 1249 1249 ! 1250 r au0 = 1026._wp !: volumic mass of reference [kg/m3]1250 rho0 = 1026._wp !: volumic mass of reference [kg/m3] 1251 1251 rcp = 3991.86795711963_wp !: heat capacity [J/K] 1252 1252 ! … … 1657 1657 WRITE(numout,*) ' ==>>> use of simplified eos: ' 1658 1658 WRITE(numout,*) ' rhd(dT=T-10,dS=S-35,Z) = [-a0*(1+lambda1/2*dT+mu1*Z)*dT ' 1659 WRITE(numout,*) ' + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / r au0'1659 WRITE(numout,*) ' + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / rho0' 1660 1660 WRITE(numout,*) ' with the following coefficients :' 1661 1661 WRITE(numout,*) ' thermal exp. coef. rn_a0 = ', rn_a0 … … 1676 1676 END SELECT 1677 1677 ! 1678 r au0_rcp = rau0 * rcp1679 r1_r au0 = 1._wp / rau01678 rho0_rcp = rho0 * rcp 1679 r1_rho0 = 1._wp / rho0 1680 1680 r1_rcp = 1._wp / rcp 1681 r1_r au0_rcp = 1._wp / rau0_rcp1681 r1_rho0_rcp = 1._wp / rho0_rcp 1682 1682 ! 1683 1683 IF(lwp) THEN … … 1694 1694 IF(lwp) WRITE(numout,*) 1695 1695 IF(lwp) WRITE(numout,*) ' Associated physical constant' 1696 IF(lwp) WRITE(numout,*) ' volumic mass of reference r au0 = ', rau0 , ' kg/m^3'1697 IF(lwp) WRITE(numout,*) ' 1. / r au0 r1_rau0 = ', r1_rau0, ' m^3/kg'1696 IF(lwp) WRITE(numout,*) ' volumic mass of reference rho0 = ', rho0 , ' kg/m^3' 1697 IF(lwp) WRITE(numout,*) ' 1. / rho0 r1_rho0 = ', r1_rho0, ' m^3/kg' 1698 1698 IF(lwp) WRITE(numout,*) ' ocean specific heat rcp = ', rcp , ' J/Kelvin' 1699 IF(lwp) WRITE(numout,*) ' r au0 * rcp rau0_rcp = ', rau0_rcp1700 IF(lwp) WRITE(numout,*) ' 1. / ( r au0 * rcp ) r1_rau0_rcp = ', r1_rau0_rcp1699 IF(lwp) WRITE(numout,*) ' rho0 * rcp rho0_rcp = ', rho0_rcp 1700 IF(lwp) WRITE(numout,*) ' 1. / ( rho0 * rcp ) r1_rho0_rcp = ', r1_rho0_rcp 1701 1701 ! 1702 1702 END SUBROUTINE eos_init -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/traadv.F90
r9598 r9939 87 87 INTEGER :: jk ! dummy loop index 88 88 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zun, zvn, zwn ! 3D workspace 89 REAL(wp), DIMENSION(:,:,: ), ALLOCATABLE :: ztrdt, ztrds89 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ztrd 90 90 !!---------------------------------------------------------------------- 91 91 ! 92 92 IF( ln_timing ) CALL timing_start('tra_adv') 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 93 ! 99 94 ! !== effective transport ==! … … 138 133 ! 139 134 IF( l_trdtra ) THEN !* Save ta and sa trends 140 ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) ) 141 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 142 ztrds(:,:,:) = tsa(:,:,:,jp_sal) 135 ALLOCATE( ztrd(jpi,jpj,jpk,jpts) ) 136 ztrd(:,:,:,:) = tsa(:,:,:,:) 143 137 ENDIF 144 138 ! … … 146 140 ! 147 141 CASE ( np_CEN ) ! Centered scheme : 2nd / 4th order 148 CALL tra_adv_cen ( kt, nit000, 'TRA', 142 CALL tra_adv_cen ( kt, nit000, 'TRA', zun, zvn, zwn , tsn, tsa, jpts, nn_cen_h, nn_cen_v ) 149 143 CASE ( np_FCT ) ! FCT scheme : 2nd / 4th order 150 CALL tra_adv_fct ( kt, nit000, 'TRA', r 2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts, nn_fct_h, nn_fct_v )144 CALL tra_adv_fct ( kt, nit000, 'TRA', rDt, zun, zvn, zwn, tsb, tsn, tsa, jpts, nn_fct_h, nn_fct_v ) 151 145 CASE ( np_MUS ) ! MUSCL 152 CALL tra_adv_mus ( kt, nit000, 'TRA', r 2dt, zun, zvn, zwn, tsb, tsa, jpts , ln_mus_ups )146 CALL tra_adv_mus ( kt, nit000, 'TRA', rDt, zun, zvn, zwn, tsb, tsa, jpts , ln_mus_ups ) 153 147 CASE ( np_UBS ) ! UBS 154 CALL tra_adv_ubs ( kt, nit000, 'TRA', r 2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts , nn_ubs_v )148 CALL tra_adv_ubs ( kt, nit000, 'TRA', rDt, zun, zvn, zwn, tsb, tsn, tsa, jpts , nn_ubs_v ) 155 149 CASE ( np_QCK ) ! QUICKEST 156 CALL tra_adv_qck ( kt, nit000, 'TRA', r 2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )150 CALL tra_adv_qck ( kt, nit000, 'TRA', rDt, zun, zvn, zwn, tsb, tsn, tsa, jpts ) 157 151 ! 158 152 END SELECT … … 160 154 IF( l_trdtra ) THEN ! save the advective trends for further diagnostics 161 155 DO jk = 1, jpkm1 162 ztrdt(:,:,jk) = tsa(:,:,jk,jp_tem) - ztrdt(:,:,jk) 163 ztrds(:,:,jk) = tsa(:,:,jk,jp_sal) - ztrds(:,:,jk) 156 ztrd(:,:,jk,:) = tsa(:,:,jk,:) - ztrd(:,:,jk,:) 164 157 END DO 165 CALL trd_tra( kt, 'TRA', jp_tem, jptra_totad, ztrd t)166 CALL trd_tra( kt, 'TRA', jp_sal, jptra_totad, ztrd s)167 DEALLOCATE( ztrd t, ztrds)158 CALL trd_tra( kt, 'TRA', jp_tem, jptra_totad, ztrd(:,:,:,jp_tem) ) 159 CALL trd_tra( kt, 'TRA', jp_sal, jptra_totad, ztrd(:,:,:,jp_sal) ) 160 DEALLOCATE( ztrd ) 168 161 ENDIF 169 162 ! ! print mean trends (used for debugging) -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/traadv_fct.F90
r9598 r9939 20 20 USE diaptr ! poleward transport diagnostics 21 21 USE diaar5 ! AR5 diagnostics 22 USE phycst , ONLY : rau0_rcp23 22 ! 24 23 USE in_out_manager ! I/O manager … … 131 130 zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) ) 132 131 zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) ) 132 !!gm faster coding ? ===>>> to be tested : 133 ! zwx(ji,jj,jk) = MAX( pun(ji,jj,jk) , 0._wp ) * ptb(ji ,jj,jk,jn) & 134 ! & + MIN( pun(ji,jj,jk) , 0._wp ) * ptb(ji+1,jj,jk,jn) 135 ! zwy(ji,jj,jk) = MAX( pvn(ji,jj,jk) , 0._wp ) * ptb(ji,jj ,jk,jn) & 136 ! & + MIN( pvn(ji,jj,jk) , 0._wp ) * ptb(ji,jj+1,jk,jn) 137 !!gm 138 133 139 END DO 134 140 END DO … … 141 147 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 142 148 zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 149 !!gm faster coding ? ===>>> to be tested : 150 ! zwx(ji,jj,jk) = MAX( pwn(ji,jj,jk) , 0._wp ) * pwn(ji,jj,jk ,jn) & 151 ! & + MIN( pwn(ji,jj,jk) , 0._wp ) * pwn(ji,jj,jk-1,jn) 152 !!gm 143 153 END DO 144 154 END DO -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/trabbc.F90
r9598 r9939 64 64 !! ocean bottom can be computed once and is added to the temperature 65 65 !! trend juste above the bottom at each time step: 66 !! ta = ta + Qsf / (r au0 rcp e3T) for k= mbkt66 !! ta = ta + Qsf / (rho0 rcp e3T) for k= mbkt 67 67 !! Where Qsf is the geothermal heat flux. 68 68 !! … … 76 76 ! 77 77 INTEGER :: ji, jj ! dummy loop indices 78 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrd t! 3D workspace78 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrd ! 3D workspace 79 79 !!---------------------------------------------------------------------- 80 80 ! … … 82 82 ! 83 83 IF( l_trdtra ) THEN ! Save the input temperature trend 84 ALLOCATE( ztrd t(jpi,jpj,jpk) )85 ztrd t(:,:,:) = tsa(:,:,:,jp_tem)84 ALLOCATE( ztrd(jpi,jpj,jpk) ) 85 ztrd(:,:,:) = tsa(:,:,:,jp_tem) 86 86 ENDIF 87 87 ! ! Add the geothermal trend on temperature … … 95 95 ! 96 96 IF( l_trdtra ) THEN ! Send the trend for diagnostics 97 ztrd t(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)98 CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbc, ztrd t)99 DEALLOCATE( ztrd t)97 ztrd(:,:,:) = tsa(:,:,:,jp_tem) - ztrd(:,:,:) 98 CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbc, ztrd ) 99 DEALLOCATE( ztrd ) 100 100 ENDIF 101 101 ! … … 157 157 ALLOCATE( qgh_trd0(jpi,jpj) ) ! allocation 158 158 ! 159 SELECT CASE ( nn_geoflx ) ! geothermal heat flux / (r auO * Cp)159 SELECT CASE ( nn_geoflx ) ! geothermal heat flux / (rhoO * Cp) 160 160 ! 161 161 CASE ( 1 ) !* constant flux 162 162 IF(lwp) WRITE(numout,*) ' ==>>> constant heat flux = ', rn_geoflx_cst 163 qgh_trd0(:,:) = r1_r au0_rcp * rn_geoflx_cst163 qgh_trd0(:,:) = r1_rho0_rcp * rn_geoflx_cst 164 164 ! 165 165 CASE ( 2 ) !* variable geothermal heat flux : read the geothermal fluxes in mW/m2 … … 178 178 179 179 CALL fld_read( nit000, 1, sf_qgh ) ! Read qgh data 180 qgh_trd0(:,:) = r1_r au0_rcp * sf_qgh(1)%fnow(:,:,1) * 1.e-3 ! conversion in W/m2180 qgh_trd0(:,:) = r1_rho0_rcp * sf_qgh(1)%fnow(:,:,1) * 1.e-3 ! conversion in W/m2 181 181 ! 182 182 CASE DEFAULT -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/trabbl.F90
r9598 r9939 103 103 INTEGER, INTENT( in ) :: kt ! ocean time-step 104 104 ! 105 REAL(wp), ALLOCATABLE, DIMENSION(:,:,: ) :: ztrdt, ztrds105 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztrd 106 106 !!---------------------------------------------------------------------- 107 107 ! … … 109 109 ! 110 110 IF( l_trdtra ) THEN !* Save the T-S input trends 111 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 112 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 113 ztrds(:,:,:) = tsa(:,:,:,jp_sal) 111 ALLOCATE( ztrd(jpi,jpj,jpk,jpts) ) 112 ztrd(:,:,:,:) = tsa(:,:,:,:) 114 113 ENDIF 115 114 … … 143 142 144 143 IF( l_trdtra ) THEN ! send the trends for further diagnostics 145 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 146 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 147 CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt ) 148 CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds ) 149 DEALLOCATE( ztrdt, ztrds ) 144 ztrd(:,:,:,:) = tsa(:,:,:,:) - ztrd(:,:,:,:) 145 CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrd(:,:,:,jp_tem) ) 146 CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrd(:,:,:,jp_sal) ) 147 DEALLOCATE( ztrd ) 150 148 ENDIF 151 149 ! -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/tradmp.F90
r9598 r9939 94 94 INTEGER :: ji, jj, jk, jn ! dummy loop indices 95 95 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts_dta 96 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ztrd ts96 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ztrd 97 97 !!---------------------------------------------------------------------- 98 98 ! … … 100 100 ! 101 101 IF( l_trdtra ) THEN !* Save ta and sa trends 102 ALLOCATE( ztrd ts(jpi,jpj,jpk,jpts) )103 ztrd ts(:,:,:,:) = tsa(:,:,:,:)102 ALLOCATE( ztrd(jpi,jpj,jpk,jpts) ) 103 ztrd(:,:,:,:) = tsa(:,:,:,:) 104 104 ENDIF 105 105 ! !== input T-S data at kt ==! … … 150 150 ! 151 151 IF( l_trdtra ) THEN ! trend diagnostic 152 ztrd ts(:,:,:,:) = tsa(:,:,:,:) - ztrdts(:,:,:,:)153 CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ztrd ts(:,:,:,jp_tem) )154 CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, ztrd ts(:,:,:,jp_sal) )155 DEALLOCATE( ztrd ts)152 ztrd(:,:,:,:) = tsa(:,:,:,:) - ztrd(:,:,:,:) 153 CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ztrd(:,:,:,jp_tem) ) 154 CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, ztrd(:,:,:,jp_sal) ) 155 DEALLOCATE( ztrd ) 156 156 ENDIF 157 157 ! ! Control print -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/traldf.F90
r9598 r9939 55 55 INTEGER, INTENT( in ) :: kt ! ocean time-step index 56 56 !! 57 REAL(wp), ALLOCATABLE, DIMENSION(:,:,: ) :: ztrdt, ztrds57 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztrd ! 4D workspace 58 58 !!---------------------------------------------------------------------- 59 59 ! … … 61 61 ! 62 62 IF( l_trdtra ) THEN !* Save ta and sa trends 63 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 64 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 65 ztrds(:,:,:) = tsa(:,:,:,jp_sal) 63 ALLOCATE( ztrd(jpi,jpj,jpk,jpts) ) 64 ztrd(:,:,:,:) = tsa(:,:,:,:) 66 65 ENDIF 67 66 ! … … 78 77 ! 79 78 IF( l_trdtra ) THEN !* save the horizontal diffusive trends for further diagnostics 80 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 81 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 82 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt ) 83 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrds ) 84 DEALLOCATE( ztrdt, ztrds ) 79 ztrd(:,:,:,:) = tsa(:,:,:,:) - ztrd(:,:,:,:) 80 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrd(:,:,:,jp_tem) ) 81 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrd(:,:,:,jp_sal) ) 82 DEALLOCATE( ztrd ) 85 83 ENDIF 86 84 ! !* print mean trends (used for debugging) -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/traldf_iso.F90
r9779 r9939 108 108 REAL(wp) :: zmsku, zahu_w, zabe1, zcof1, zcoef3 ! local scalars 109 109 REAL(wp) :: zmskv, zahv_w, zabe2, zcof2, zcoef4 ! - - 110 REAL(wp) :: zcoef0, ze3w_2, zsign , z2dt, z1_2dt! - -110 REAL(wp) :: zcoef0, ze3w_2, zsign ! - - 111 111 REAL(wp), DIMENSION(jpi,jpj) :: zdkt, zdk1t, z2d 112 112 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdit, zdjt, zftu, zftv, ztfw … … 127 127 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 128 128 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 129 !130 ! ! set time step size (Euler/Leapfrog)131 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2dt = rdt ! at nit000 (Euler)132 ELSE ; z2dt = 2.* rdt ! (Leapfrog)133 ENDIF134 z1_2dt = 1._wp / z2dt135 129 ! 136 130 IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0) … … 191 185 DO ji = 1, fs_jpim1 192 186 ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 193 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 )194 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 195 189 END DO 196 190 END DO -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/traldf_triad.F90
r9598 r9939 85 85 INTEGER :: ip,jp,kp ! dummy loop indices 86 86 INTEGER :: ierr ! local integer 87 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 88 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 89 REAL(wp) :: zcoef0, ze3w_2, zsign , z2dt, z1_2dt! - -87 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars 88 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 89 REAL(wp) :: zcoef0, ze3w_2, zsign ! - - 90 90 ! 91 91 REAL(wp) :: zslope_skew, zslope_iso, zslope2, zbu, zbv … … 110 110 l_hst = .FALSE. 111 111 l_ptr = .FALSE. 112 IF( cdtype == 'TRA' .AND. ln_diaptr ) l_ptr = .TRUE. 113 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 114 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 115 ! 116 ! ! set time step size (Euler/Leapfrog) 117 IF( neuler == 0 .AND. kt == kit000 ) THEN ; z2dt = rdt ! at nit000 (Euler) 118 ELSE ; z2dt = 2.* rdt ! (Leapfrog) 112 IF( cdtype == 'TRA' ) THEN 113 IF ( ln_diaptr ) l_ptr = .TRUE. 114 IF ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 115 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) l_hst = .TRUE. 119 116 ENDIF 120 z1_2dt = 1._wp / z2dt121 117 ! 122 118 IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0) … … 202 198 DO ji = 1, fs_jpim1 203 199 ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 204 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 )205 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt200 zcoef0 = rDt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 201 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_Dt 206 202 END DO 207 203 END DO -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/tramle.F90
r9598 r9939 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 … … 115 115 zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 116 116 zmld(ji,jj) = zmld(ji,jj) + zc 117 zbm (ji,jj) = zbm (ji,jj) + zc * (r au0 - rhop(ji,jj,jk) ) * r1_rau0117 zbm (ji,jj) = zbm (ji,jj) + zc * (rho0 - rhop(ji,jj,jk) ) * r1_rho0 118 118 zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 119 119 END DO … … 302 302 IF( ln_mle ) THEN ! MLE initialisation 303 303 ! 304 rb_c = grav * rn_rho_c_mle / rau0! Mixed Layer buoyancy criteria304 rb_c = grav * rn_rho_c_mle / rho0 ! Mixed Layer buoyancy criteria 305 305 IF(lwp) WRITE(numout,*) 306 306 IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 ' -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/tranpc.F90
r9598 r9939 65 65 LOGICAL :: l_bottom_reached, l_column_treated 66 66 REAL(wp) :: zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z 67 REAL(wp) :: zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw , z1_r2dt67 REAL(wp) :: zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw 68 68 REAL(wp), PARAMETER :: zn2_zero = 1.e-14_wp ! acceptance criteria for neutrality (N2==0) 69 69 REAL(wp), DIMENSION( jpk ) :: zvn2 ! vertical profile of N2 at 1 given point... … … 71 71 REAL(wp), DIMENSION(jpi,jpj,jpk ) :: zn2 ! N^2 72 72 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zab ! alpha and beta 73 REAL(wp), ALLOCATABLE, DIMENSION(:,:,: ) :: ztrdt, ztrds ! 3D workspace73 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztrd ! 4D workspace 74 74 ! 75 75 LOGICAL, PARAMETER :: l_LB_debug = .FALSE. ! set to true if you want to follow what is … … 82 82 IF( MOD( kt, nn_npc ) == 0 ) THEN 83 83 ! 84 IF( l_trdtra ) THEN !* Save initial after fields 85 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 86 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 87 ztrds(:,:,:) = tsa(:,:,:,jp_sal) 84 IF( l_trdtra ) THEN !* Save input after fields 85 ALLOCATE( ztrd(jpi,jpj,jpk,jpts) ) 86 ztrd(:,:,:,:) = tsa(:,:,:,:) 88 87 ENDIF 89 88 ! … … 301 300 ! 302 301 IF( l_trdtra ) THEN ! send the Non penetrative mixing trends for diagnostic 303 z1_r2dt = 1._wp / (2._wp * rdt) 304 ztrdt(:,:,:) = ( tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ) * z1_r2dt 305 ztrds(:,:,:) = ( tsa(:,:,:,jp_sal) - ztrds(:,:,:) ) * z1_r2dt 306 CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt ) 307 CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds ) 308 DEALLOCATE( ztrdt, ztrds ) 302 ztrd(:,:,:,:) = ( tsa(:,:,:,:) - ztrd(:,:,:,:) ) * r1_Dt 303 CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrd(:,:,:,jp_tem) ) 304 CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrd(:,:,:,jp_sal) ) 305 DEALLOCATE( ztrd ) 309 306 ENDIF 310 307 ! -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/tranxt.F90
r9598 r9939 90 90 INTEGER :: ji, jj, jk, jn ! dummy loop indices 91 91 REAL(wp) :: zfact ! local scalars 92 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, ztrds92 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ztrd ! 4D workspace 93 93 !!---------------------------------------------------------------------- 94 94 ! … … 111 111 IF( ln_bdy ) CALL bdy_tra( kt ) ! BDY open boundaries 112 112 113 ! set time step size (Euler/Leapfrog) 114 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! at nit000 (Euler) 115 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2._wp* rdt ! at nit000 or nit000+1 (Leapfrog) 116 ENDIF 117 118 ! trends computation initialisation 119 IF( l_trdtra ) THEN 120 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 121 ztrdt(:,:,jpk) = 0._wp 122 ztrds(:,:,jpk) = 0._wp 113 IF( l_trdtra ) THEN ! trends computation initialisation 114 ALLOCATE( ztrd(jpi,jpj,jpk,jpts) ) 115 ztrd(:,:,jpk,:) = 0._wp 123 116 IF( ln_traldf_iso ) THEN ! diagnose the "pure" Kz diffusive trend 124 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrd t)125 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrd s)117 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrd(:,:,:,jp_tem) ) 118 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrd(:,:,:,jp_sal) ) 126 119 ENDIF 127 120 ! total trend for the non-time-filtered variables. 128 zfact = 1.0 / r dt121 zfact = 1.0 / rn_Dt 129 122 ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from tsn terms 130 DO jk = 1, jpkm1 131 ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) / e3t_n(:,:,jk) - tsn(:,:,jk,jp_tem)) * zfact 132 ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) / e3t_n(:,:,jk) - tsn(:,:,jk,jp_sal)) * zfact 133 END DO 134 CALL trd_tra( kt, 'TRA', jp_tem, jptra_tot, ztrdt ) 135 CALL trd_tra( kt, 'TRA', jp_sal, jptra_tot, ztrds ) 136 IF( ln_linssh ) THEN ! linear sea surface height only 137 ! Store now fields before applying the Asselin filter 138 ! in order to calculate Asselin filter trend later. 139 ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 140 ztrds(:,:,:) = tsn(:,:,:,jp_sal) 141 ENDIF 142 ENDIF 143 144 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step (only swap) 123 DO jn = 1, jpts 124 DO jk = 1, jpkm1 125 ztrd(:,:,jk,jn) = ( tsa(:,:,jk,jn)*e3t_a(:,:,jk) / e3t_n(:,:,jk) - tsn(:,:,jk,jn)) * zfact 126 END DO 127 END DO 128 CALL trd_tra( kt, 'TRA', jp_tem, jptra_tot, ztrd(:,:,:,jp_tem) ) 129 CALL trd_tra( kt, 'TRA', jp_sal, jptra_tot, ztrd(:,:,:,jp_sal) ) 130 IF( ln_linssh ) THEN ! linear sea surface height only Store now fields before applying 131 ! ! the Asselin filter in order to calculate Asselin filter trend later. 132 ztrd(:,:,:,:) = tsn(:,:,:,:) 133 ENDIF 134 ENDIF 135 136 IF( l_1st_euler ) THEN ! Euler time-stepping at first time-step (only swap) 145 137 DO jn = 1, jpts 146 138 DO jk = 1, jpkm1 … … 150 142 IF (l_trdtra .AND. .NOT. ln_linssh ) THEN ! Zero Asselin filter contribution must be explicitly written out since for vvl 151 143 ! ! Asselin filter is output by tra_nxt_vvl that is not called on this time step 152 ztrdt(:,:,:) = 0._wp 153 ztrds(:,:,:) = 0._wp 154 CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt ) 155 CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds ) 144 ztrd(:,:,:,:) = 0._wp 145 CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrd(:,:,:,jp_tem) ) 146 CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrd(:,:,:,jp_sal) ) 156 147 END IF 157 148 ! 158 ELSE ! Leap-Frog + Asselin filter time stepping 159 ! 160 IF( ln_linssh ) THEN ; CALL tra_nxt_fix( kt, nit000, 'TRA', tsb, tsn, tsa, jpts ) ! linear free surface 161 ELSE ; CALL tra_nxt_vvl( kt, nit000, rdt, 'TRA', tsb, tsn, tsa, & 162 & sbc_tsc, sbc_tsc_b, jpts ) ! non-linear free surface 163 ENDIF 164 ! 165 CALL lbc_lnk_multi( tsb(:,:,:,jp_tem), 'T', 1., tsb(:,:,:,jp_sal), 'T', 1., & 166 & tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., & 167 & tsa(:,:,:,jp_tem), 'T', 1., tsa(:,:,:,jp_sal), 'T', 1. ) 149 ELSE ! Leap-Frog + Asselin filter time stepping 150 ! 151 IF( ln_linssh ) THEN ; CALL tra_nxt_fix( kt, nit000, 'TRA', tsb, tsn, tsa, jpts ) ! linear free surface 152 ELSE ; CALL tra_nxt_vvl( kt, nit000, rn_Dt,'TRA', tsb, tsn, tsa, & 153 & sbc_tsc, sbc_tsc_b, jpts ) ! non-linear free surface 154 ENDIF 155 ! 156 CALL lbc_lnk_multi( tsb, 'T', 1., tsn, 'T', 1., tsa, 'T', 1. ) 168 157 ! 169 158 ENDIF 170 159 ! 171 160 IF( l_trdtra .AND. ln_linssh ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt 172 zfact = 1._wp / r2dt173 161 DO jk = 1, jpkm1 174 ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact 175 ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact 176 END DO 177 CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt ) 178 CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds ) 162 ztrd(:,:,jk,:) = ( tsb(:,:,jk,:) - ztrd(:,:,jk,:) ) * r1_Dt 163 END DO 164 CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrd(:,:,:,jp_tem) ) 165 CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrd(:,:,:,jp_sal) ) 179 166 END IF 180 IF( l_trdtra ) DEALLOCATE( ztrd t , ztrds)167 IF( l_trdtra ) DEALLOCATE( ztrd ) 181 168 ! 182 169 ! ! control print … … 227 214 ztd = pta(ji,jj,jk,jn) - 2._wp * ztn + ptb(ji,jj,jk,jn) ! time laplacian on tracers 228 215 ! 229 ptb(ji,jj,jk,jn) = ztn + atfp * ztd! ptb <-- filtered ptn216 ptb(ji,jj,jk,jn) = ztn + rn_atfp * ztd ! ptb <-- filtered ptn 230 217 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 231 218 END DO … … 238 225 239 226 240 SUBROUTINE tra_nxt_vvl( kt, kit000, p 2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )227 SUBROUTINE tra_nxt_vvl( kt, kit000, pdt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt ) 241 228 !!---------------------------------------------------------------------- 242 229 !! *** ROUTINE tra_nxt_vvl *** … … 247 234 !! ** Method : - Apply a thickness weighted Asselin time filter on now fields. 248 235 !! - swap tracer fields to prepare the next time_step. 249 !! tb = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )250 !! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] )236 !! tb = ( e3t_n*tn + rn_atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 237 !! /( e3t_n + rn_atfp*[ e3t_b - 2 e3t_n + e3t_a ] ) 251 238 !! tn = ta 252 239 !! … … 255 242 INTEGER , INTENT(in ) :: kt ! ocean time-step index 256 243 INTEGER , INTENT(in ) :: kit000 ! first time step index 257 REAL(wp) , INTENT(in ) :: p 2dt! time-step244 REAL(wp) , INTENT(in ) :: pdt ! time-step 258 245 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 259 246 INTEGER , INTENT(in ) :: kjpt ! number of tracers … … 289 276 IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) ) THEN 290 277 ALLOCATE( ztrd_atf(jpi,jpj,jpk,kjpt) ) 291 ztrd_atf(:,:,:,:) = 0.0_wp 292 ENDIF 293 zfact = 1._wp / r2dt 294 zfact1 = atfp * p2dt 295 zfact2 = zfact1 * r1_rau0 296 DO jn = 1, kjpt 278 ztrd_atf(:,:,:,:) = 0._wp 279 ENDIF 280 ! 281 zfact = r1_Dt 282 zfact1 = rn_atfp * pdt 283 zfact2 = zfact1 * r1_rho0 284 DO jn = 1, kjpt 297 285 DO jk = 1, jpkm1 298 286 DO jj = 2, jpjm1 … … 309 297 ztc_d = ztc_a - 2. * ztc_n + ztc_b 310 298 ! 311 ze3t_f = ze3t_n + atfp * ze3t_d312 ztc_f = ztc_n + atfp * ztc_d299 ze3t_f = ze3t_n + rn_atfp * ze3t_d 300 ztc_f = ztc_n + rn_atfp * ztc_d 313 301 ! 314 302 IF( jk == mikt(ji,jj) ) THEN ! first level -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/traqsr.F90
r9598 r9939 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 … … 112 112 REAL(wp) :: zlogc, zlogc2, zlogc3 113 113 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zekb, zekg, zekr 114 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrd t114 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrd 115 115 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zetot, zchl3d 116 116 !!---------------------------------------------------------------------- … … 125 125 ! 126 126 IF( l_trdtra ) THEN ! trends diagnostic: save the input temperature trend 127 ALLOCATE( ztrd t(jpi,jpj,jpk) )128 ztrd t(:,:,:) = tsa(:,:,:,jp_tem)127 ALLOCATE( ztrd(jpi,jpj,jpk) ) 128 ztrd(:,:,:) = tsa(:,:,:,jp_tem) 129 129 ENDIF 130 130 ! … … 133 133 ! !-----------------------------------! 134 134 IF( kt == nit000 ) THEN !== 1st time step ==! 135 !!gm case neuler not taken into account.... 136 IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN ! read in restart 135 IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 .AND. .NOT.l_1st_euler ) THEN ! read in restart 137 136 IF(lwp) WRITE(numout,*) ' nit000-1 qsr tracer content forcing field read in the restart file' 138 137 z1_2 = 0.5_wp … … 154 153 ! 155 154 DO jk = 1, nksr 156 qsr_hc(:,:,jk) = r1_r au0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )155 qsr_hc(:,:,jk) = r1_rho0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) ) 157 156 END DO 158 157 ! … … 234 233 DO jj = 2, jpjm1 235 234 DO ji = fs_2, fs_jpim1 236 qsr_hc(ji,jj,jk) = r1_r au0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) )235 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) ) 237 236 END DO 238 237 END DO … … 243 242 CASE( np_2BD ) !== 2-bands fluxes ==! 244 243 ! 245 zz0 = rn_abs * r1_r au0_rcp ! surface equi-partition in 2-bands246 zz1 = ( 1. - rn_abs ) * r1_r au0_rcp244 zz0 = rn_abs * r1_rho0_rcp ! surface equi-partition in 2-bands 245 zz1 = ( 1. - rn_abs ) * r1_rho0_rcp 247 246 DO jk = 1, nksr ! solar heat absorbed at T-point in the top 400m 248 247 DO jj = 2, jpjm1 … … 270 269 DO jj = 2, jpjm1 271 270 DO ji = fs_2, fs_jpim1 ! vector opt. 272 IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_r au0_rcp * qsr(ji,jj) )271 IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) ) 273 272 ELSE ; fraqsr_1lev(ji,jj) = 1._wp 274 273 ENDIF … … 281 280 zetot(:,:,nksr+1:jpk) = 0._wp ! below ~400m set to zero 282 281 DO jk = nksr, 1, -1 283 zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * r au0_rcp282 zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp 284 283 END DO 285 284 CALL iom_put( 'qsr3d', zetot ) ! 3D distribution of shortwave Radiation … … 295 294 ! 296 295 IF( l_trdtra ) THEN ! qsr tracers trends saved for diagnostics 297 ztrd t(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)298 CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrd t)299 DEALLOCATE( ztrd t)296 ztrd(:,:,:) = tsa(:,:,:,jp_tem) - ztrd(:,:,:) 297 CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrd ) 298 DEALLOCATE( ztrd ) 300 299 ENDIF 301 300 ! ! print mean trends (used for debugging) -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/trasbc.F90
r9598 r9939 78 78 INTEGER :: ikt, ikb ! local integers 79 79 REAL(wp) :: zfact, z1_e3t, zdep, ztim ! local scalar 80 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, ztrds80 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ztrd ! 4D workspace 81 81 !!---------------------------------------------------------------------- 82 82 ! … … 89 89 ENDIF 90 90 ! 91 IF( l_trdtra ) THEN !* Save ta and sa trends 92 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 93 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 94 ztrds(:,:,:) = tsa(:,:,:,jp_sal) 91 IF( l_trdtra ) THEN !* Save input tsa trends 92 ALLOCATE( ztrd(jpi,jpj,jpk,jpts) ) 93 ztrd(:,:,:,:) = tsa(:,:,:,:) 95 94 ENDIF 96 95 ! … … 98 97 IF( .NOT.ln_traqsr ) THEN ! no solar radiation penetration 99 98 qns(:,:) = qns(:,:) + qsr(:,:) ! total heat flux in qns 100 qsr(:,:) = 0._wp 99 qsr(:,:) = 0._wp ! qsr set to zero 101 100 ENDIF 102 101 … … 127 126 IF ( ll_wd ) THEN ! If near WAD point limit the flux for now 128 127 IF ( sshn(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 129 sbc_tsc(ji,jj,jp_tem) = r1_r au0_rcp * qns(ji,jj) ! non solar heat flux128 sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj) ! non solar heat flux 130 129 ELSE IF ( sshn(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 131 sbc_tsc(ji,jj,jp_tem) = r1_r au0_rcp * qns(ji,jj) &130 sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj) & 132 131 & * tanh ( 5._wp * ( ( sshn(ji,jj) + ht_0(ji,jj) - rn_wdmin1 ) * r_rn_wdmin1 ) ) 133 132 ELSE … … 135 134 ENDIF 136 135 ELSE 137 sbc_tsc(ji,jj,jp_tem) = r1_r au0_rcp * qns(ji,jj) ! non solar heat flux136 sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj) ! non solar heat flux 138 137 ENDIF 139 138 140 sbc_tsc(ji,jj,jp_sal) = r1_r au0 * sfx(ji,jj) ! salt flux due to freezing/melting139 sbc_tsc(ji,jj,jp_sal) = r1_rho0 * sfx(ji,jj) ! salt flux due to freezing/melting 141 140 END DO 142 141 END DO … … 144 143 DO jj = 2, jpj !==>> add concentration/dilution effect due to constant volume cell 145 144 DO ji = fs_2, fs_jpim1 ! vector opt. 146 sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_r au0 * emp(ji,jj) * tsn(ji,jj,1,jp_tem)147 sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_r au0 * emp(ji,jj) * tsn(ji,jj,1,jp_sal)145 sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rho0 * emp(ji,jj) * tsn(ji,jj,1,jp_tem) 146 sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rho0 * emp(ji,jj) * tsn(ji,jj,1,jp_sal) 148 147 END DO 149 148 END DO !==>> output c./d. term … … 272 271 273 272 IF( l_trdtra ) THEN ! save the horizontal diffusive trends for further diagnostics 274 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 275 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 276 CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt ) 277 CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds ) 278 DEALLOCATE( ztrdt , ztrds ) 273 ztrd(:,:,:,:) = tsa(:,:,:,:) - ztrd(:,:,:,:) 274 CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrd(:,:,:,jp_tem) ) 275 CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrd(:,:,:,jp_sal) ) 276 DEALLOCATE( ztrd ) 279 277 ENDIF 280 278 ! -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/trazdf.F90
r9598 r9939 52 52 INTEGER, INTENT(in) :: kt ! ocean time-step index 53 53 ! 54 INTEGER :: jk ! Dummy loop indices55 REAL(wp), DIMENSION(:,:,: ), ALLOCATABLE :: ztrdt, ztrds ! 3D workspace54 INTEGER :: jk, jts ! Dummy loop indices 55 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ztrd ! 4D workspace 56 56 !!--------------------------------------------------------------------- 57 57 ! … … 64 64 ENDIF 65 65 ! 66 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! at nit000, = rdt (restarting with Euler time stepping) 67 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2. * rdt ! otherwise, = 2 rdt (leapfrog) 66 IF( l_trdtra ) THEN !* Save input tsa trend 67 ALLOCATE( ztrd(jpi,jpj,jpk,jpts) ) 68 ztrd(:,:,:,:) = tsa(:,:,:,:) 68 69 ENDIF 69 70 ! 70 IF( l_trdtra ) THEN !* Save ta and sa trends71 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )72 ztrdt(:,:,:) = tsa(:,:,:,jp_tem)73 ztrds(:,:,:) = tsa(:,:,:,jp_sal)74 ENDIF75 !76 71 ! !* compute lateral mixing trend and add it to the general trend 77 CALL tra_zdf_imp( kt, nit000, 'TRA', r 2dt, tsb, tsa, jpts )72 CALL tra_zdf_imp( kt, nit000, 'TRA', rDt, tsb, tsa, jpts ) 78 73 79 74 !!gm WHY here ! and I don't like that ! … … 85 80 86 81 IF( l_trdtra ) THEN ! save the vertical diffusive trends for further diagnostics 87 DO j k = 1, jpkm188 ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_tem)*e3t_b(:,:,jk) ) &89 & / (e3t_n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk)90 ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_sal)*e3t_b(:,:,jk) ) &91 & / (e3t_n(:,:,jk)*r2dt) ) - ztrds(:,:,jk)82 DO jts = 1, jpts 83 DO jk = 1, jpkm1 84 ztrd(:,:,jk,jts) = ( ( tsa(:,:,jk,jts)*e3t_a(:,:,jk) - tsb(:,:,jk,jts)*e3t_b(:,:,jk) ) / (e3t_n(:,:,jk)*rDt) ) & 85 & - ztrd(:,:,jk,jts) 86 END DO 92 87 END DO 93 88 !!gm this should be moved in trdtra.F90 and done on all trends 94 CALL lbc_lnk _multi( ztrdt, 'T', 1. , ztrds, 'T', 1. )89 CALL lbc_lnk( ztrd, 'T', 1. ) 95 90 !!gm 96 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrd t)97 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrd s)98 DEALLOCATE( ztrd t , ztrds)91 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrd(:,:,:,jp_tem) ) 92 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrd(:,:,:,jp_sal) ) 93 DEALLOCATE( ztrd ) 99 94 ENDIF 100 95 ! ! print mean trends (used for debugging) … … 180 175 DO jj = 2, jpjm1 181 176 DO ji = fs_2, fs_jpim1 ! vector opt. 182 !!gm BUG I think, use e3w_a instead of e3w_n, not sure of that183 177 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w_n(ji,jj,jk ) 184 178 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)
Note: See TracChangeset
for help on using the changeset viewer.
