- Timestamp:
- 2020-04-08T21:37:59+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3
- Files:
-
- 14 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3
- 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/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/eosbn2.F90
r12622 r12724 192 192 !! *** ROUTINE eos_insitu *** 193 193 !! 194 !! ** Purpose : Compute the in situ density (ratio rho/r au0) from194 !! ** Purpose : Compute the in situ density (ratio rho/rho0) from 195 195 !! potential temperature and salinity using an equation of state 196 196 !! selected in the nameos namelist 197 197 !! 198 !! ** Method : prd(t,s,z) = ( rho(t,s,z) - r au0 ) / rau0198 !! ** Method : prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0 199 199 !! with prd in situ density anomaly no units 200 200 !! t TEOS10: CT or EOS80: PT Celsius … … 202 202 !! z depth meters 203 203 !! rho in situ density kg/m^3 204 !! r au0 reference density kg/m^3204 !! rho0 reference density kg/m^3 205 205 !! 206 206 !! ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z). … … 211 211 !! 212 212 !! ln_seos : simplified equation of state 213 !! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / r au0213 !! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0 214 214 !! linear case function of T only: rn_alpha<>0, other coefficients = 0 215 215 !! 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_3D … … 284 284 & - rn_nu * zt * zs 285 285 ! 286 prd(ji,jj,jk) = zn * r1_r au0 * ztm ! density anomaly (masked)286 prd(ji,jj,jk) = zn * r1_rho0 * ztm ! density anomaly (masked) 287 287 END_3D 288 288 ! … … 300 300 !! *** ROUTINE eos_insitu_pot *** 301 301 !! 302 !! ** Purpose : Compute the in situ density (ratio rho/r au0) and the302 !! ** Purpose : Compute the in situ density (ratio rho/rho0) and the 303 303 !! potential volumic mass (Kg/m3) from potential temperature and 304 304 !! salinity fields using an equation of state selected in the … … 380 380 prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp) ! potential density referenced at the surface 381 381 ! 382 prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_r au0 - 1._wp ) ! density anomaly (masked)382 prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_rho0 - 1._wp ) ! density anomaly (masked) 383 383 END DO 384 384 prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos … … 420 420 prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface 421 421 ! 422 prd(ji,jj,jk) = ( zn * r1_r au0 - 1._wp ) * ztm ! density anomaly (masked)422 prd(ji,jj,jk) = ( zn * r1_rho0 - 1._wp ) * ztm ! density anomaly (masked) 423 423 END_3D 424 424 ENDIF … … 435 435 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs & 436 436 & - rn_nu * zt * zs 437 prhop(ji,jj,jk) = ( r au0 + zn ) * ztm437 prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 438 438 ! ! density anomaly (masked) 439 439 zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 440 prd(ji,jj,jk) = zn * r1_r au0 * ztm440 prd(ji,jj,jk) = zn * r1_rho0 * ztm 441 441 ! 442 442 END_3D … … 455 455 !! *** ROUTINE eos_insitu_2d *** 456 456 !! 457 !! ** Purpose : Compute the in situ density (ratio rho/r au0) from457 !! ** Purpose : Compute the in situ density (ratio rho/rho0) from 458 458 !! potential temperature and salinity using an equation of state 459 459 !! selected in the nameos namelist. * 2D field case … … 509 509 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 510 510 ! 511 prd(ji,jj) = zn * r1_r au0 - 1._wp ! unmasked in situ density anomaly511 prd(ji,jj) = zn * r1_rho0 - 1._wp ! unmasked in situ density anomaly 512 512 ! 513 513 END_2D … … 525 525 & - rn_nu * zt * zs 526 526 ! 527 prd(ji,jj) = zn * r1_r au0 ! unmasked in situ density anomaly527 prd(ji,jj) = zn * r1_rho0 ! unmasked in situ density anomaly 528 528 ! 529 529 END_2D … … 589 589 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 590 590 ! 591 pab(ji,jj,jk,jp_tem) = zn * r1_r au0 * ztm591 pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm 592 592 ! 593 593 ! beta … … 610 610 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 611 611 ! 612 pab(ji,jj,jk,jp_sal) = zn / zs * r1_r au0 * ztm612 pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm 613 613 ! 614 614 END_3D … … 623 623 ! 624 624 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 625 pab(ji,jj,jk,jp_tem) = zn * r1_r au0 * ztm ! alpha625 pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm ! alpha 626 626 ! 627 627 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 628 pab(ji,jj,jk,jp_sal) = zn * r1_r au0 * ztm ! beta628 pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm ! beta 629 629 ! 630 630 END_3D … … 695 695 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 696 696 ! 697 pab(ji,jj,jp_tem) = zn * r1_r au0697 pab(ji,jj,jp_tem) = zn * r1_rho0 698 698 ! 699 699 ! beta … … 716 716 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 717 717 ! 718 pab(ji,jj,jp_sal) = zn / zs * r1_r au0718 pab(ji,jj,jp_sal) = zn / zs * r1_rho0 719 719 ! 720 720 ! … … 730 730 ! 731 731 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 732 pab(ji,jj,jp_tem) = zn * r1_r au0 ! alpha732 pab(ji,jj,jp_tem) = zn * r1_rho0 ! alpha 733 733 ! 734 734 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 735 pab(ji,jj,jp_sal) = zn * r1_r au0 ! beta735 pab(ji,jj,jp_sal) = zn * r1_rho0 ! beta 736 736 ! 737 737 END_2D … … 800 800 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 801 801 ! 802 pab(jp_tem) = zn * r1_r au0802 pab(jp_tem) = zn * r1_rho0 803 803 ! 804 804 ! beta … … 821 821 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 822 822 ! 823 pab(jp_sal) = zn / zs * r1_r au0823 pab(jp_sal) = zn / zs * r1_rho0 824 824 ! 825 825 ! … … 832 832 ! 833 833 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 834 pab(jp_tem) = zn * r1_r au0 ! alpha834 pab(jp_tem) = zn * r1_rho0 ! alpha 835 835 ! 836 836 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 837 pab(jp_sal) = zn * r1_r au0 ! beta837 pab(jp_sal) = zn * r1_rho0 ! beta 838 838 ! 839 839 CASE DEFAULT … … 1053 1053 !! ** Method : PE is defined analytically as the vertical 1054 1054 !! primitive of EOS times -g integrated between 0 and z>0. 1055 !! pen is the nonlinear bsq-PE anomaly: pen = ( PE - r au0 gz ) / rau0 gz - rd1055 !! pen is the nonlinear bsq-PE anomaly: pen = ( PE - rho0 gz ) / rho0 gz - rd 1056 1056 !! = 1/z * /int_0^z rd dz - rd 1057 1057 !! where rd is the density anomaly (see eos_rhd function) 1058 1058 !! ab_pe are partial derivatives of PE anomaly with respect to T and S: 1059 !! ab_pe(1) = - 1/(r au0 gz) * dPE/dT + drd/dT = - d(pen)/dT1060 !! ab_pe(2) = 1/(r au0 gz) * dPE/dS + drd/dS = d(pen)/dS1059 !! ab_pe(1) = - 1/(rho0 gz) * dPE/dT + drd/dT = - d(pen)/dT 1060 !! ab_pe(2) = 1/(rho0 gz) * dPE/dS + drd/dS = d(pen)/dS 1061 1061 !! 1062 1062 !! ** Action : - pen : PE anomaly given at T-points … … 1104 1104 zn = ( zn2 * zh + zn1 ) * zh + zn0 1105 1105 ! 1106 ppen(ji,jj,jk) = zn * zh * r1_r au0 * ztm1106 ppen(ji,jj,jk) = zn * zh * r1_rho0 * ztm 1107 1107 ! 1108 1108 ! alphaPE non-linear anomaly … … 1119 1119 zn = ( zn2 * zh + zn1 ) * zh + zn0 1120 1120 ! 1121 pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_r au0 * ztm1121 pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm 1122 1122 ! 1123 1123 ! betaPE non-linear anomaly … … 1134 1134 zn = ( zn2 * zh + zn1 ) * zh + zn0 1135 1135 ! 1136 pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_r au0 * ztm1136 pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm 1137 1137 ! 1138 1138 END_3D … … 1145 1145 zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point 1146 1146 ztm = tmask(ji,jj,jk) ! tmask 1147 zn = 0.5_wp * zh * r1_r au0 * ztm1147 zn = 0.5_wp * zh * r1_rho0 * ztm 1148 1148 ! ! Potential Energy 1149 1149 ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn … … 1187 1187 IF(lwm) WRITE( numond, nameos ) 1188 1188 ! 1189 r au0 = 1026._wp !: volumic mass of reference [kg/m3]1189 rho0 = 1026._wp !: volumic mass of reference [kg/m3] 1190 1190 rcp = 3991.86795711963_wp !: heat capacity [J/K] 1191 1191 ! … … 1599 1599 WRITE(numout,*) ' ==>>> use of simplified eos: ' 1600 1600 WRITE(numout,*) ' rhd(dT=T-10,dS=S-35,Z) = [-a0*(1+lambda1/2*dT+mu1*Z)*dT ' 1601 WRITE(numout,*) ' + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / r au0'1601 WRITE(numout,*) ' + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / rho0' 1602 1602 WRITE(numout,*) ' with the following coefficients :' 1603 1603 WRITE(numout,*) ' thermal exp. coef. rn_a0 = ', rn_a0 … … 1618 1618 END SELECT 1619 1619 ! 1620 r au0_rcp = rau0 * rcp1621 r1_r au0 = 1._wp / rau01620 rho0_rcp = rho0 * rcp 1621 r1_rho0 = 1._wp / rho0 1622 1622 r1_rcp = 1._wp / rcp 1623 r1_r au0_rcp = 1._wp / rau0_rcp1623 r1_rho0_rcp = 1._wp / rho0_rcp 1624 1624 ! 1625 1625 IF(lwp) THEN … … 1636 1636 IF(lwp) WRITE(numout,*) 1637 1637 IF(lwp) WRITE(numout,*) ' Associated physical constant' 1638 IF(lwp) WRITE(numout,*) ' volumic mass of reference r au0 = ', rau0 , ' kg/m^3'1639 IF(lwp) WRITE(numout,*) ' 1. / r au0 r1_rau0 = ', r1_rau0, ' m^3/kg'1638 IF(lwp) WRITE(numout,*) ' volumic mass of reference rho0 = ', rho0 , ' kg/m^3' 1639 IF(lwp) WRITE(numout,*) ' 1. / rho0 r1_rho0 = ', r1_rho0, ' m^3/kg' 1640 1640 IF(lwp) WRITE(numout,*) ' ocean specific heat rcp = ', rcp , ' J/Kelvin' 1641 IF(lwp) WRITE(numout,*) ' r au0 * rcp rau0_rcp = ', rau0_rcp1642 IF(lwp) WRITE(numout,*) ' 1. / ( r au0 * rcp ) r1_rau0_rcp = ', r1_rau0_rcp1641 IF(lwp) WRITE(numout,*) ' rho0 * rcp rho0_rcp = ', rho0_rcp 1642 IF(lwp) WRITE(numout,*) ' 1. / ( rho0 * rcp ) r1_rho0_rcp = ', r1_rho0_rcp 1643 1643 ! 1644 1644 END SUBROUTINE eos_init -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv.F90
r12624 r12724 93 93 IF( ln_timing ) CALL timing_start('tra_adv') 94 94 ! 95 ! ! set time step96 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! at nit000 (Euler)97 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2._wp * rdt ! at nit000 or nit000+1 (Leapfrog)98 ENDIF99 !100 95 ! !== effective transport ==! 101 96 zuu(:,:,jpk) = 0._wp … … 153 148 CALL tra_adv_cen ( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v ) 154 149 CASE ( np_FCT ) ! FCT scheme : 2nd / 4th order 155 CALL tra_adv_fct ( kt, nit000, 'TRA', r 2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v )150 CALL tra_adv_fct ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v ) 156 151 CASE ( np_MUS ) ! MUSCL 157 CALL tra_adv_mus ( kt, nit000, 'TRA', r 2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )152 CALL tra_adv_mus ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 158 153 CASE ( np_UBS ) ! UBS 159 CALL tra_adv_ubs ( kt, nit000, 'TRA', r 2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v )154 CALL tra_adv_ubs ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v ) 160 155 CASE ( np_QCK ) ! QUICKEST 161 CALL tra_adv_qck ( kt, nit000, 'TRA', r 2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs )156 CALL tra_adv_qck ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs ) 162 157 ! 163 158 END SELECT -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv_fct.F90
r12590 r12724 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/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatf.F90
r12680 r12724 114 114 IF( ln_bdy ) CALL bdy_tra( kt, Kbb, pts, Kaa ) ! BDY open boundaries 115 115 116 ! set time step size (Euler/Leapfrog)117 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! at nit000 (Euler)118 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2._wp* rdt ! at nit000 or nit000+1 (Leapfrog)119 ENDIF120 121 116 ! trends computation initialisation 122 117 IF( l_trdtra ) THEN … … 129 124 ENDIF 130 125 ! total trend for the non-time-filtered variables. 131 zfact = 1.0 / r dt126 zfact = 1.0 / rn_Dt 132 127 ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from pts(Kmm) terms 133 128 DO jk = 1, jpkm1 … … 145 140 ENDIF 146 141 147 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping142 IF( l_1st_euler ) THEN ! Euler time-stepping 148 143 ! 149 144 IF (l_trdtra .AND. .NOT. ln_linssh ) THEN ! Zero Asselin filter contribution must be explicitly written out since for vvl … … 157 152 ELSE ! Leap-Frog + Asselin filter time stepping 158 153 ! 159 IF( ln_linssh ) THEN ; CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000, 'TRA', pts, jpts ) ! linear free surface160 ELSE ; CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nit000, r dt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts ) ! non-linear free surface154 IF( ln_linssh ) THEN ; CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000, 'TRA', pts, jpts ) ! linear free surface 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 161 156 ENDIF 162 157 ! … … 167 162 ENDIF 168 163 ! 169 IF( l_trdtra .AND. ln_linssh ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt 170 zfact = 1._wp / r2dt 164 IF( l_trdtra .AND. ln_linssh ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt 171 165 DO jk = 1, jpkm1 172 ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kmm) - ztrdt(:,:,jk) ) * zfact173 ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kmm) - ztrds(:,:,jk) ) * zfact166 ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kmm) - ztrdt(:,:,jk) ) * r1_Dt 167 ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kmm) - ztrds(:,:,jk) ) * r1_Dt 174 168 END DO 175 169 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_atf, ztrdt ) … … 220 214 ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb) ! time laplacian on tracers 221 215 ! 222 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 223 217 END_3D 224 218 ! … … 235 229 !! 236 230 !! ** Method : - Apply a thickness weighted Asselin time filter on now fields. 237 !! pt(Kmm) = ( e3t_Kmm*pt(Kmm) + atfp*[ e3t_Kbb*pt(Kbb) - 2 e3t_Kmm*pt(Kmm) + e3t_Kaa*pt(Kaa) ] )238 !! /( 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_Kaa*pt(Kaa) ] ) 232 !! /( e3t_Kmm + rn_atfp*[ e3t_Kbb - 2 e3t_Kmm + e3t_Kaa ] ) 239 233 !! 240 234 !! ** Action : - pt(Kmm) ready for the next time step … … 278 272 ENDIF 279 273 zfact = 1._wp / p2dt 280 zfact1 = atfp * p2dt281 zfact2 = zfact1 * r1_r au0274 zfact1 = rn_atfp * p2dt 275 zfact2 = zfact1 * r1_rho0 282 276 DO jn = 1, kjpt 283 277 DO_3D_00_00( 1, jpkm1 ) … … 293 287 ztc_d = ztc_a - 2. * ztc_n + ztc_b 294 288 ! 295 ze3t_f = ze3t_n + atfp * ze3t_d296 ztc_f = ztc_n + atfp * ztc_d289 ze3t_f = ze3t_n + rn_atfp * ze3t_d 290 ztc_f = ztc_n + rn_atfp * ztc_d 297 291 ! 298 292 ! Add asselin correction on scale factors: -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatfQCO.F90
r12624 r12724 114 114 ! IF( ln_bdy ) CALL bdy_tra( kt, Kbb, pts, Kaa ) ! BDY open boundaries 115 115 116 ! set time step size (Euler/Leapfrog)117 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! at nit000 (Euler)118 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2._wp* rdt ! at nit000 or nit000+1 (Leapfrog)119 ENDIF120 121 116 ! trends computation initialisation 122 117 IF( l_trdtra ) THEN … … 129 124 ENDIF 130 125 ! total trend for the non-time-filtered variables. 131 zfact = 1.0 / r dt126 zfact = 1.0 / rn_Dt 132 127 ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from pts(Kmm) terms 133 128 DO jk = 1, jpkm1 … … 149 144 ENDIF 150 145 151 IF( neuler == 0 .AND. kt == nit000) THEN ! Euler time-stepping146 IF( l_1st_euler ) THEN ! Euler time-stepping 152 147 ! 153 148 IF (l_trdtra .AND. .NOT. ln_linssh ) THEN ! Zero Asselin filter contribution must be explicitly written out since for vvl … … 161 156 ELSE ! Leap-Frog + Asselin filter time stepping 162 157 ! 163 IF ( ln_linssh ) THEN ; CALL tra_atf_fix_lf( kt, Kbb, Kmm, Kaa, nit000, 'TRA', pts, jpts ) ! linear free surface164 ELSE ; CALL tra_atf_qco_lf( kt, Kbb, Kmm, Kaa, nit000, r dt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts ) ! non-linear free surface158 IF ( ln_linssh ) THEN ; CALL tra_atf_fix_lf( kt, Kbb, Kmm, Kaa, nit000, 'TRA', pts, jpts ) ! linear free surface 159 ELSE ; CALL tra_atf_qco_lf( kt, Kbb, Kmm, Kaa, nit000, rn_Dt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts ) ! non-linear free surface 165 160 ENDIF 166 161 ! … … 172 167 ! 173 168 IF( l_trdtra .AND. ln_linssh ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt 174 zfact = 1._wp / r2dt175 169 DO jk = 1, jpkm1 176 ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kmm) - ztrdt(:,:,jk) ) * zfact177 ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kmm) - ztrds(:,:,jk) ) * zfact170 ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kmm) - ztrdt(:,:,jk) ) * r1_Dt 171 ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kmm) - ztrds(:,:,jk) ) * r1_Dt 178 172 END DO 179 173 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_atf, ztrdt ) … … 224 218 ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb) ! time laplacian on tracers 225 219 ! 226 pt(ji,jj,jk,jn,Kmm) = ztn + atfp * ztd ! pt <-- filtered pt220 pt(ji,jj,jk,jn,Kmm) = ztn + rn_atfp * ztd ! pt <-- filtered pt 227 221 END_3D 228 222 ! … … 239 233 !! 240 234 !! ** Method : - Apply a thickness weighted Asselin time filter on now fields. 241 !! pt(Kmm) = ( e3t_m*pt(Kmm) + atfp*[ e3t_b*pt(Kbb) - 2 e3t_m*pt(Kmm) + e3t_a*pt(Kaa) ] )242 !! /( e3t_m + atfp*[ e3t_b - 2 e3t_m + e3t_a ] )235 !! pt(Kmm) = ( e3t_m*pt(Kmm) + rn_atfp*[ e3t_b*pt(Kbb) - 2 e3t_m*pt(Kmm) + e3t_a*pt(Kaa) ] ) 236 !! /( e3t_m + rn_atfp*[ e3t_b - 2 e3t_m + e3t_a ] ) 243 237 !! 244 238 !! ** Action : - pt(Kmm) ready for the next time step … … 282 276 ENDIF 283 277 zfact = 1._wp / p2dt 284 zfact1 = atfp * p2dt285 zfact2 = zfact1 * r1_r au0278 zfact1 = rn_atfp * p2dt 279 zfact2 = zfact1 * r1_rho0 286 280 DO jn = 1, kjpt 287 281 DO_3D_00_00( 1, jpkm1 ) … … 297 291 ztc_d = ztc_a - 2. * ztc_n + ztc_b 298 292 ! 299 ztc_f = ztc_n + atfp * ztc_d293 ztc_f = ztc_n + rn_atfp * ztc_d 300 294 ! 301 295 ! Asselin correction on scale factors is done via ssh in r3t_f -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/trabbc.F90
r12590 r12724 67 67 !! ocean bottom can be computed once and is added to the temperature 68 68 !! trend juste above the bottom at each time step: 69 !! ta = ta + Qsf / (r au0 rcp e3T) for k= mbkt69 !! ta = ta + Qsf / (rho0 rcp e3T) for k= mbkt 70 70 !! Where Qsf is the geothermal heat flux. 71 71 !! … … 104 104 ENDIF 105 105 ! 106 CALL iom_put ( "hfgeou" , r au0_rcp * qgh_trd0(:,:) )106 CALL iom_put ( "hfgeou" , rho0_rcp * qgh_trd0(:,:) ) 107 107 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbc - Ta: ', mask1=tmask, clinfo3='tra-ta' ) 108 108 ! … … 164 164 CASE ( 1 ) !* constant flux 165 165 IF(lwp) WRITE(numout,*) ' ==>>> constant heat flux = ', rn_geoflx_cst 166 qgh_trd0(:,:) = r1_r au0_rcp * rn_geoflx_cst166 qgh_trd0(:,:) = r1_rho0_rcp * rn_geoflx_cst 167 167 ! 168 168 CASE ( 2 ) !* variable geothermal heat flux : read the geothermal fluxes in mW/m2 … … 181 181 182 182 CALL fld_read( nit000, 1, sf_qgh ) ! Read qgh data 183 qgh_trd0(:,:) = r1_r au0_rcp * sf_qgh(1)%fnow(:,:,1) * 1.e-3 ! conversion in W/m2183 qgh_trd0(:,:) = r1_rho0_rcp * sf_qgh(1)%fnow(:,:,1) * 1.e-3 ! conversion in W/m2 184 184 ! 185 185 CASE DEFAULT -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traldf_iso.F90
r12622 r12724 110 110 REAL(wp) :: zmsku, zahu_w, zabe1, zcof1, zcoef3 ! local scalars 111 111 REAL(wp) :: zmskv, zahv_w, zabe2, zcof2, zcoef4 ! - - 112 REAL(wp) :: zcoef0, ze3w_2, zsign , z2dt, z1_2dt! - -112 REAL(wp) :: zcoef0, ze3w_2, zsign ! - - 113 113 REAL(wp), DIMENSION(jpi,jpj) :: zdkt, zdk1t, z2d 114 114 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdit, zdjt, zftu, zftv, ztfw … … 130 130 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 131 131 ! 132 ! ! set time step size (Euler/Leapfrog)133 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2dt = rdt ! at nit000 (Euler)134 ELSE ; z2dt = 2.* rdt ! (Leapfrog)135 ENDIF136 z1_2dt = 1._wp / z2dt137 132 ! 138 133 IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0) … … 182 177 DO_3D_10_10( 2, jpkm1 ) 183 178 ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 184 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 )185 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt179 zcoef0 = rDt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 180 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_Dt 186 181 END_3D 187 182 ENDIF -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traldf_triad.F90
r12622 r12724 87 87 INTEGER :: ip,jp,kp ! dummy loop indices 88 88 INTEGER :: ierr ! local integer 89 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 90 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 91 REAL(wp) :: zcoef0, ze3w_2, zsign , z2dt, z1_2dt! - -89 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars 90 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 91 REAL(wp) :: zcoef0, ze3w_2, zsign ! - - 92 92 ! 93 93 REAL(wp) :: zslope_skew, zslope_iso, zslope2, zbu, zbv … … 112 112 l_hst = .FALSE. 113 113 l_ptr = .FALSE. 114 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) ) l_ptr = .TRUE. 115 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 116 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 117 ! 118 ! ! set time step size (Euler/Leapfrog) 119 IF( neuler == 0 .AND. kt == kit000 ) THEN ; z2dt = rdt ! at nit000 (Euler) 120 ELSE ; z2dt = 2.* rdt ! (Leapfrog) 114 IF( cdtype == 'TRA' ) THEN 115 IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf') ) l_ptr = .TRUE. 116 IF( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 117 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) l_hst = .TRUE. 121 118 ENDIF 122 z1_2dt = 1._wp / z2dt123 119 ! 124 120 IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0) … … 193 189 DO_3D_10_10( 2, jpkm1 ) 194 190 ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 195 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 )196 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt191 zcoef0 = rDt * ( 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 * r1_Dt 197 193 END_3D 198 194 ENDIF -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/tramle.F90
r12590 r12724 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 … … 113 113 zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 114 114 zmld(ji,jj) = zmld(ji,jj) + zc 115 zbm (ji,jj) = zbm (ji,jj) + zc * (r au0 - rhop(ji,jj,jk) ) * r1_rau0115 zbm (ji,jj) = zbm (ji,jj) + zc * (rho0 - rhop(ji,jj,jk) ) * r1_rho0 116 116 zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 117 117 END_3D … … 274 274 IF( ln_mle ) THEN ! MLE initialisation 275 275 ! 276 rb_c = grav * rn_rho_c_mle /r au0 ! Mixed Layer buoyancy criteria276 rb_c = grav * rn_rho_c_mle /rho0 ! Mixed Layer buoyancy criteria 277 277 IF(lwp) WRITE(numout,*) 278 278 IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 ' -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/tranpc.F90
r12590 r12724 68 68 LOGICAL :: l_bottom_reached, l_column_treated 69 69 REAL(wp) :: zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z 70 REAL(wp) :: zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r 2dt70 REAL(wp) :: zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_rDt 71 71 REAL(wp), PARAMETER :: zn2_zero = 1.e-14_wp ! acceptance criteria for neutrality (N2==0) 72 72 REAL(wp), DIMENSION( jpk ) :: zvn2 ! vertical profile of N2 at 1 given point... … … 302 302 ! 303 303 IF( l_trdtra ) THEN ! send the Non penetrative mixing trends for diagnostic 304 z1_r 2dt = 1._wp / (2._wp * rdt)305 ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * z1_r 2dt306 ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * z1_r 2dt304 z1_rDt = 1._wp / (2._wp * rn_Dt) 305 ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * z1_rDt 306 ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * z1_rDt 307 307 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_npc, ztrdt ) 308 308 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_npc, ztrds ) -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traqsr.F90
r12590 r12724 88 88 !! I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 89 89 !! The temperature trend associated with the solar radiation penetration 90 !! is given by : zta = 1/e3t dk[ I ] / (r au0*Cp)90 !! is given by : zta = 1/e3t dk[ I ] / (rho0*Cp) 91 91 !! At the bottom, boudary condition for the radiation is no flux : 92 92 !! all heat which has not been absorbed in the above levels is put … … 136 136 ! !-----------------------------------! 137 137 IF( kt == nit000 ) THEN !== 1st time step ==! 138 !!gm case neuler not taken into account.... 139 IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN ! read in restart 138 IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 .AND. .NOT.l_1st_euler ) THEN ! read in restart 140 139 IF(lwp) WRITE(numout,*) ' nit000-1 qsr tracer content forcing field read in the restart file' 141 140 z1_2 = 0.5_wp … … 157 156 ! 158 157 DO jk = 1, nksr 159 qsr_hc(:,:,jk) = r1_r au0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )158 qsr_hc(:,:,jk) = r1_rho0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) ) 160 159 END DO 161 160 ! … … 229 228 ! 230 229 DO_3D_00_00( 1, nksr ) 231 qsr_hc(ji,jj,jk) = r1_r au0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) )230 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) ) 232 231 END_3D 233 232 ! … … 236 235 CASE( np_2BD ) !== 2-bands fluxes ==! 237 236 ! 238 zz0 = rn_abs * r1_r au0_rcp ! surface equi-partition in 2-bands239 zz1 = ( 1. - rn_abs ) * r1_r au0_rcp237 zz0 = rn_abs * r1_rho0_rcp ! surface equi-partition in 2-bands 238 zz1 = ( 1. - rn_abs ) * r1_rho0_rcp 240 239 DO_3D_00_00( 1, nksr ) 241 240 zc0 = zz0 * EXP( -gdepw(ji,jj,jk ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk ,Kmm)*xsi1r ) … … 255 254 ! sea-ice: store the 1st ocean level attenuation coefficient 256 255 DO_2D_00_00 257 IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_r au0_rcp * qsr(ji,jj) )256 IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) ) 258 257 ELSE ; fraqsr_1lev(ji,jj) = 1._wp 259 258 ENDIF … … 265 264 zetot(:,:,nksr+1:jpk) = 0._wp ! below ~400m set to zero 266 265 DO jk = nksr, 1, -1 267 zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * r au0_rcp266 zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp 268 267 END DO 269 268 CALL iom_put( 'qsr3d', zetot ) ! 3D distribution of shortwave Radiation -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/trasbc.F90
r12590 r12724 125 125 ! !== Now sbc tracer content fields ==! 126 126 DO_2D_01_00 127 sbc_tsc(ji,jj,jp_tem) = r1_r au0_rcp * qns(ji,jj) ! non solar heat flux128 sbc_tsc(ji,jj,jp_sal) = r1_r au0 * sfx(ji,jj) ! salt flux due to freezing/melting127 sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj) ! non solar heat flux 128 sbc_tsc(ji,jj,jp_sal) = r1_rho0 * sfx(ji,jj) ! salt flux due to freezing/melting 129 129 END_2D 130 130 IF( ln_linssh ) THEN !* linear free surface 131 131 DO_2D_01_00 132 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)133 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)132 sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm) 133 sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm) 134 134 END_2D 135 135 IF( iom_use('emp_x_sst') ) CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) ) -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/trazdf.F90
r12590 r12724 67 67 ENDIF 68 68 ! 69 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! at nit000, = rdt (restarting with Euler time stepping)70 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2. * rdt ! otherwise, = 2 rdt (leapfrog)71 ENDIF72 !73 69 IF( l_trdtra ) THEN !* Save ta and sa trends 74 70 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) … … 78 74 ! 79 75 ! !* compute lateral mixing trend and add it to the general trend 80 CALL tra_zdf_imp( kt, nit000, 'TRA', r 2dt, Kbb, Kmm, Krhs, pts, Kaa, jpts )76 CALL tra_zdf_imp( kt, nit000, 'TRA', rDt, Kbb, Kmm, Krhs, pts, Kaa, jpts ) 81 77 82 78 !!gm WHY here ! and I don't like that ! … … 89 85 IF( l_trdtra ) THEN ! save the vertical diffusive trends for further diagnostics 90 86 DO jk = 1, jpkm1 91 ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) & 92 & - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) & 93 & / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrdt(:,:,jk) 94 ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) & 95 & - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) & 96 & / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrds(:,:,jk) 87 ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) & 88 & - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) & 89 & / ( e3t(:,:,jk,Kmm)*rDt ) ) & 90 & - ztrdt(:,:,jk) 91 ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) & 92 & - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) & 93 & / ( e3t(:,:,jk,Kmm)*rDt ) ) & 94 & - ztrds(:,:,jk) 97 95 END DO 98 96 !!gm this should be moved in trdtra.F90 and done on all trends
Note: See TracChangeset
for help on using the changeset viewer.