Changeset 15066
- Timestamp:
- 2021-07-01T13:17:50+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE
- Files:
-
- 8 edited
-
C1D/step_c1d.F90 (modified) (2 diffs)
-
DOM/domqco.F90 (modified) (2 diffs)
-
DYN/dynspg_ts.F90 (modified) (5 diffs)
-
TRA/eosbn2.F90 (modified) (14 diffs)
-
TRA/tranpc.F90 (modified) (1 diff)
-
ZDF/zdfosm.F90 (modified) (1 diff)
-
step.F90 (modified) (1 diff)
-
stpmlf.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/C1D/step_c1d.F90
r14644 r15066 72 72 CALL eos_rab( CASTWP(ts(:,:,:,:,Nbb)), rab_b, Nnn ) ! before local thermal/haline expension ratio at T-points 73 73 CALL eos_rab( CASTWP(ts(:,:,:,:,Nnn)), rab_n, Nnn ) ! now local thermal/haline expension ratio at T-points 74 CALL bn2( CASTWP(ts(:,:,:,:,Nbb)), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency75 CALL bn2( CASTWP(ts(:,:,:,:,Nnn)), rab_n, rn2 , Nnn ) ! now Brunt-Vaisala frequency74 CALL bn2( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency 75 CALL bn2( ts(:,:,:,:,Nnn), rab_n, rn2 , Nnn ) ! now Brunt-Vaisala frequency 76 76 77 77 ! VERTICAL PHYSICS … … 108 108 IF( ln_zdfosm ) CALL tra_osm( kstp, Nnn , ts, Nrhs ) ! OSMOSIS non-local tracer fluxes 109 109 CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa ) ! vertical mixing 110 CALL eos( CAST EWP(ts(:,:,:,:,Nnn)), rhd, rhop, gdept_0(:,:,:) ) ! now potential density for zdfmxl110 CALL eos( CASTWP(ts(:,:,:,:,Nnn)), rhd, rhop, gdept_0(:,:,:) ) ! now potential density for zdfmxl 111 111 IF( ln_zdfnpc ) CALL tra_npc( kstp, Nnn, Nrhs, ts, Naa ) ! applied non penetrative convective adjustment on (t,s) 112 112 CALL tra_atf( kstp, Nbb, Nnn, Naa, ts ) ! time filtering of "now" tracer arrays -
NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DOM/domqco.F90
r14986 r15066 118 118 ! ! Horizontal interpolation of e3t 119 119 #if defined key_RK3 120 CALL dom_qco_r3c( CASTWP(ssh(:,:,Kbb)), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb), r3f(:,:) )121 CALL dom_qco_r3c( CASTWP(ssh(:,:,Kmm)), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm) )120 CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb), r3f(:,:) ) 121 CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm) ) 122 122 #else 123 CALL dom_qco_r3c( CASTWP(ssh(:,:,Kbb)), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb) )124 CALL dom_qco_r3c( CASTWP(ssh(:,:,Kmm)), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) )123 CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb) ) 124 CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) ) 125 125 #endif 126 126 ! dom_qco_r3c defines over [nn_hls, nn_hls-1, nn_hls, nn_hls-1] … … 142 142 !! - compute the ratio ssh/h_0 at t-,u-,v-pts, (f-pt optional) 143 143 !!---------------------------------------------------------------------- 144 REAL( wp), DIMENSION(:,:) , INTENT(in ) :: pssh ! sea surface height [m]144 REAL(dp), DIMENSION(:,:) , INTENT(in ) :: pssh ! sea surface height [m] 145 145 REAL(wp), DIMENSION(:,:) , INTENT( out) :: pr3t, pr3u, pr3v ! ssh/h0 ratio at t-, u-, v-,points [-] 146 146 REAL(wp), DIMENSION(:,:), OPTIONAL, INTENT( out) :: pr3f ! ssh/h0 ratio at f-point [-] -
NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DYN/dynspg_ts.F90
r14986 r15066 87 87 88 88 89 INTERFACE dyn_drg_init90 MODULE PROCEDURE dyn_drg_init_wp, dyn_drg_init_mixed91 END INTERFACE92 93 89 !! * Substitutions 94 90 # include "do_loop_substitute.h90" … … 288 284 END_2D 289 285 ELSE !* remove baroclinic drag AND provide the barotropic drag coefficients 290 CALL dyn_drg_init( Kbb, Kmm, CASTWP(puu), CASTWP(pvv), puu_b, pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v )286 CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b, pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v ) 291 287 ENDIF 292 288 ! … … 1347 1343 1348 1344 1349 SUBROUTINE dyn_drg_init _wp( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )1345 SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v ) 1350 1346 !!---------------------------------------------------------------------- 1351 1347 !! *** ROUTINE dyn_drg_init *** … … 1358 1354 !!---------------------------------------------------------------------- 1359 1355 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 1360 REAL( wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in ) :: puu, pvv ! ocean velocities and RHS of momentum equation1356 REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in ) :: puu, pvv ! ocean velocities and RHS of momentum equation 1361 1357 REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(in ) :: puu_b, pvv_b ! barotropic velocities at main time levels 1362 1358 REAL(wp), DIMENSION(jpi,jpj) , INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS … … 1451 1447 ENDIF 1452 1448 ! 1453 END SUBROUTINE dyn_drg_init_wp 1454 1455 SUBROUTINE dyn_drg_init_mixed( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v ) 1456 !!---------------------------------------------------------------------- 1457 !! *** ROUTINE dyn_drg_init *** 1458 !! 1459 !! ** Purpose : - add the baroclinic top/bottom drag contribution to 1460 !! the baroclinic part of the barotropic RHS 1461 !! - compute the barotropic drag coefficients 1462 !! 1463 !! ** Method : computation done over the INNER domain only 1464 !!---------------------------------------------------------------------- 1465 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 1466 REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in ) :: puu, pvv ! ocean velocities and RHS of momentum equation 1467 REAL(sp), DIMENSION(jpi,jpj,jpt) , INTENT(in ) :: puu_b, pvv_b ! barotropic velocities at main time levels 1468 REAL(sp), DIMENSION(jpi,jpj) , INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS 1469 REAL(sp), DIMENSION(jpi,jpj) , INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients 1470 ! 1471 INTEGER :: ji, jj ! dummy loop indices 1472 INTEGER :: ikbu, ikbv, iktu, iktv 1473 REAL(sp) :: zztmp 1474 REAL(sp), DIMENSION(jpi,jpj) :: zu_i, zv_i 1475 !!---------------------------------------------------------------------- 1476 ! 1477 ! !== Set the barotropic drag coef. ==! 1478 ! 1479 IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! top+bottom friction (ocean cavities) 1480 1481 DO_2D( 0, 0, 0, 0 ) 1482 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 1483 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 1484 END_2D 1485 ELSE ! bottom friction only 1486 DO_2D( 0, 0, 0, 0 ) 1487 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 1488 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 1489 END_2D 1490 ENDIF 1491 ! 1492 ! !== BOTTOM stress contribution from baroclinic velocities ==! 1493 ! 1494 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities 1495 1496 DO_2D( 0, 0, 0, 0 ) 1497 ikbu = mbku(ji,jj) 1498 ikbv = mbkv(ji,jj) 1499 zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm) 1500 zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm) 1501 END_2D 1502 ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities 1503 1504 DO_2D( 0, 0, 0, 0 ) 1505 ikbu = mbku(ji,jj) 1506 ikbv = mbkv(ji,jj) 1507 zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb) 1508 zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb) 1509 END_2D 1510 ENDIF 1511 ! 1512 IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please ! 1513 zztmp = -1._wp / rDt_e 1514 DO_2D( 0, 0, 0, 0 ) 1515 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1516 & r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) 1517 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( & 1518 & r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) 1519 END_2D 1520 ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 1521 1522 DO_2D( 0, 0, 0, 0 ) 1523 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj) 1524 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj) 1525 END_2D 1526 END IF 1527 ! 1528 ! !== TOP stress contribution from baroclinic velocities ==! (no W/D case) 1529 ! 1530 IF( ln_isfcav.OR.ln_drgice_imp ) THEN 1531 ! 1532 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity 1533 1534 DO_2D( 0, 0, 0, 0 ) 1535 iktu = miku(ji,jj) 1536 iktv = mikv(ji,jj) 1537 zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm) 1538 zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm) 1539 END_2D 1540 ELSE ! CENTRED integration: use BEFORE top baroclinic velocity 1541 1542 DO_2D( 0, 0, 0, 0 ) 1543 iktu = miku(ji,jj) 1544 iktv = mikv(ji,jj) 1545 zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb) 1546 zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb) 1547 END_2D 1548 ENDIF 1549 ! 1550 ! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 1551 1552 DO_2D( 0, 0, 0, 0 ) 1553 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) 1554 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj) 1555 END_2D 1556 ! 1557 ENDIF 1558 ! 1559 END SUBROUTINE dyn_drg_init_mixed 1449 END SUBROUTINE dyn_drg_init 1560 1450 1561 1451 SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in -
NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/eosbn2.F90
r14986 r15066 56 56 ! !! * Interface 57 57 INTERFACE eos 58 MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d, eos_insitu_pot_2d 58 MODULE PROCEDURE eos_insitu_wp, eos_insitu_pot_wp, eos_insitu_2d, eos_insitu_pot_2d 59 MODULE PROCEDURE eos_insitu_mixed, eos_insitu_pot_mixed 59 60 END INTERFACE 60 61 ! 61 62 INTERFACE eos_rab 62 MODULE PROCEDURE rab_3d, rab_2d, rab_0d 63 MODULE PROCEDURE rab_3d_wp, rab_2d, rab_0d 64 MODULE PROCEDURE rab_3d_mixed 63 65 END INTERFACE 64 66 ! … … 190 192 CONTAINS 191 193 192 SUBROUTINE eos_insitu ( pts, prd, pdep )194 SUBROUTINE eos_insitu_wp( pts, prd, pdep ) 193 195 !! 194 196 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] … … 197 199 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdep ! depth [m] 198 200 !! 199 CALL eos_insitu_t( pts, is_tile(pts), prd, is_tile(prd), pdep, is_tile(pdep) ) 200 END SUBROUTINE eos_insitu 201 202 SUBROUTINE eos_insitu_t( pts, ktts, prd, ktrd, pdep, ktdep ) 203 !!---------------------------------------------------------------------- 204 !! *** ROUTINE eos_insitu *** 201 CALL eos_insitu_t_wp( pts, is_tile(pts), prd, is_tile(prd), pdep, is_tile(pdep) ) 202 END SUBROUTINE eos_insitu_wp 203 204 SUBROUTINE eos_insitu_mixed( pts, prd, pdep ) 205 !! 206 REAL(dp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 207 ! ! 2 : salinity [psu] 208 REAL(sp), DIMENSION(:,:,:) , INTENT( out) :: prd ! in situ density [-] 209 REAL(sp), DIMENSION(:,:,:) , INTENT(in ) :: pdep ! depth [m] 210 !! 211 CALL eos_insitu_t_mixed( pts, is_tile(pts), prd, is_tile(prd), pdep, is_tile(pdep) ) 212 END SUBROUTINE 213 214 SUBROUTINE eos_insitu_t_wp( pts, ktts, prd, ktrd, pdep, ktdep ) 215 !!---------------------------------------------------------------------- 216 !! *** ROUTINE 205 217 !! 206 218 !! ** Purpose : Compute the in situ density (ratio rho/rho0) from … … 306 318 IF( ln_timing ) CALL timing_stop('eos-insitu') 307 319 ! 308 END SUBROUTINE eos_insitu_t 309 310 311 SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep ) 320 END SUBROUTINE eos_insitu_t_wp 321 SUBROUTINE eos_insitu_t_mixed( pts, ktts, prd, ktrd, pdep, ktdep ) 322 !!---------------------------------------------------------------------- 323 !! *** ROUTINE eos_insitu *** 324 !! 325 !! ** Purpose : Compute the in situ density (ratio rho/rho0) from 326 !! potential temperature and salinity using an equation of state 327 !! selected in the nameos namelist 328 !! 329 !! ** Method : prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0 330 !! with prd in situ density anomaly no units 331 !! t TEOS10: CT or EOS80: PT Celsius 332 !! s TEOS10: SA or EOS80: SP TEOS10: g/kg or EOS80: psu 333 !! z depth meters 334 !! rho in situ density kg/m^3 335 !! rho0 reference density kg/m^3 336 !! 337 !! ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z). 338 !! Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celsius, sa=35.5 g/kg 339 !! 340 !! ln_eos80 : polynomial EOS-80 equation of state is used for rho(t,s,z). 341 !! Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celsius, sp=35.5 psu 342 !! 343 !! ln_seos : simplified equation of state 344 !! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0 345 !! linear case function of T only: rn_alpha<>0, other coefficients = 0 346 !! linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0 347 !! Vallis like equation: use default values of coefficients 348 !! 349 !! ** Action : compute prd , the in situ density (no units) 350 !! 351 !! References : Roquet et al, Ocean Modelling, in preparation (2014) 352 !! Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 353 !! TEOS-10 Manual, 2010 354 !!---------------------------------------------------------------------- 355 INTEGER , INTENT(in ) :: ktts, ktrd, ktdep 356 REAL(dp), DIMENSION(A2D_T(ktts) ,JPK,JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 357 ! ! 2 : salinity [psu] 358 REAL(sp), DIMENSION(A2D_T(ktrd) ,JPK ), INTENT( out) :: prd ! in situ density [-] 359 REAL(sp), DIMENSION(A2D_T(ktdep),JPK ), INTENT(in ) :: pdep ! depth [m] 360 ! 361 INTEGER :: ji, jj, jk ! dummy loop indices 362 REAL(sp) :: zt , zh , zs , ztm ! local scalars 363 REAL(sp) :: zn , zn0, zn1, zn2, zn3 ! - - 364 !!---------------------------------------------------------------------- 365 ! 366 IF( ln_timing ) CALL timing_start('eos-insitu') 367 ! 368 SELECT CASE( neos ) 369 ! 370 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 371 ! 372 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 373 ! 374 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 375 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 376 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 377 ztm = tmask(ji,jj,jk) ! tmask 378 ! 379 zn3 = EOS013*zt & 380 & + EOS103*zs+EOS003 381 ! 382 zn2 = (EOS022*zt & 383 & + EOS112*zs+EOS012)*zt & 384 & + (EOS202*zs+EOS102)*zs+EOS002 385 ! 386 zn1 = (((EOS041*zt & 387 & + EOS131*zs+EOS031)*zt & 388 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 389 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 390 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 391 ! 392 zn0 = (((((EOS060*zt & 393 & + EOS150*zs+EOS050)*zt & 394 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 395 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 396 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 397 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 398 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 399 ! 400 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 401 ! 402 prd(ji,jj,jk) = ( zn * r1_rho0 - 1._wp ) * ztm ! density anomaly (masked) 403 ! 404 END_3D 405 ! 406 CASE( np_seos ) !== simplified EOS ==! 407 ! 408 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 409 zt = pts (ji,jj,jk,jp_tem) - 10._wp 410 zs = pts (ji,jj,jk,jp_sal) - 35._wp 411 zh = pdep (ji,jj,jk) 412 ztm = tmask(ji,jj,jk) 413 ! 414 zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt & 415 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs & 416 & - rn_nu * zt * zs 417 ! 418 prd(ji,jj,jk) = zn * r1_rho0 * ztm ! density anomaly (masked) 419 END_3D 420 ! 421 END SELECT 422 ! 423 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=REAL(prd, wp), clinfo1=' eos-insitu : ', kdim=jpk ) 424 ! 425 IF( ln_timing ) CALL timing_stop('eos-insitu') 426 ! 427 END SUBROUTINE eos_insitu_t_mixed 428 429 430 SUBROUTINE eos_insitu_pot_wp( pts, prd, prhop, pdep ) 312 431 !! 313 432 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] … … 317 436 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdep ! depth [m] 318 437 !! 319 CALL eos_insitu_pot_t( pts, is_tile(pts), prd, is_tile(prd), prhop, is_tile(prhop), pdep, is_tile(pdep) ) 320 END SUBROUTINE eos_insitu_pot 321 322 323 SUBROUTINE eos_insitu_pot_t( pts, ktts, prd, ktrd, prhop, ktrhop, pdep, ktdep ) 438 CALL eos_insitu_pot_t_wp( pts, is_tile(pts), prd, is_tile(prd), prhop, is_tile(prhop), pdep, is_tile(pdep) ) 439 END SUBROUTINE eos_insitu_pot_wp 440 441 SUBROUTINE eos_insitu_pot_mixed( pts, prd, prhop, pdep ) 442 !! 443 REAL(dp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 444 ! ! 2 : salinity [psu] 445 REAL(sp), DIMENSION(:,:,:) , INTENT( out) :: prd ! in situ density [-] 446 REAL(dp), DIMENSION(:,:,:) , INTENT( out) :: prhop ! potential density (surface referenced) 447 REAL(sp), DIMENSION(:,:,:) , INTENT(in ) :: pdep ! depth [m] 448 !! 449 CALL eos_insitu_pot_t_mixed( pts, is_tile(pts), prd, is_tile(prd), prhop, is_tile(prhop), pdep, is_tile(pdep) ) 450 END SUBROUTINE eos_insitu_pot_mixed 451 452 453 SUBROUTINE eos_insitu_pot_t_wp( pts, ktts, prd, ktrd, prhop, ktrhop, pdep, ktdep ) 324 454 !!---------------------------------------------------------------------- 325 455 !! *** ROUTINE eos_insitu_pot *** … … 475 605 IF( ln_timing ) CALL timing_stop('eos-pot') 476 606 ! 477 END SUBROUTINE eos_insitu_pot_t 607 END SUBROUTINE eos_insitu_pot_t_wp 608 609 SUBROUTINE eos_insitu_pot_t_mixed( pts, ktts, prd, ktrd, prhop, ktrhop, pdep, ktdep ) 610 !!---------------------------------------------------------------------- 611 !! *** ROUTINE eos_insitu_pot *** 612 !! 613 !! ** Purpose : Compute the in situ density (ratio rho/rho0) and the 614 !! potential volumic mass (Kg/m3) from potential temperature and 615 !! salinity fields using an equation of state selected in the 616 !! namelist. 617 !! 618 !! ** Action : - prd , the in situ density (no units) 619 !! - prhop, the potential volumic mass (Kg/m3) 620 !! 621 !!---------------------------------------------------------------------- 622 INTEGER , INTENT(in ) :: ktts, ktrd, ktrhop, ktdep 623 REAL(dp), DIMENSION(A2D_T(ktts) ,JPK,JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 624 ! ! 2 : salinity [psu] 625 REAL(sp), DIMENSION(A2D_T(ktrd) ,JPK ), INTENT( out) :: prd ! in situ density [-] 626 REAL(dp), DIMENSION(A2D_T(ktrhop),JPK ), INTENT( out) :: prhop ! potential density (surface referenced) 627 REAL(sp), DIMENSION(A2D_T(ktdep) ,JPK ), INTENT(in ) :: pdep ! depth [m] 628 ! 629 INTEGER :: ji, jj, jk, jsmp ! dummy loop indices 630 INTEGER :: jdof 631 REAL(wp) :: zt , zh , zstemp, zs , ztm ! local scalars 632 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 633 REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign ! local vectors 634 !!---------------------------------------------------------------------- 635 ! 636 IF( ln_timing ) CALL timing_start('eos-pot') 637 ! 638 SELECT CASE ( neos ) 639 ! 640 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 641 ! 642 ! Stochastic equation of state 643 IF ( ln_sto_eos ) THEN 644 ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 645 ALLOCATE(zn_sto(1:2*nn_sto_eos)) 646 ALLOCATE(zsign(1:2*nn_sto_eos)) 647 DO jsmp = 1, 2*nn_sto_eos, 2 648 zsign(jsmp) = 1._wp 649 zsign(jsmp+1) = -1._wp 650 END DO 651 ! 652 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 653 ! 654 ! compute density (2*nn_sto_eos) times: 655 ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 656 ! (2) for t-dt, s-ds (with the opposite fluctuation) 657 DO jsmp = 1, nn_sto_eos*2 658 jdof = (jsmp + 1) / 2 659 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 660 zt = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0 ! temperature 661 zstemp = pts (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 662 zs = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 ) ! square root salinity 663 ztm = tmask(ji,jj,jk) ! tmask 664 ! 665 zn3 = EOS013*zt & 666 & + EOS103*zs+EOS003 667 ! 668 zn2 = (EOS022*zt & 669 & + EOS112*zs+EOS012)*zt & 670 & + (EOS202*zs+EOS102)*zs+EOS002 671 ! 672 zn1 = (((EOS041*zt & 673 & + EOS131*zs+EOS031)*zt & 674 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 675 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 676 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 677 ! 678 zn0_sto(jsmp) = (((((EOS060*zt & 679 & + EOS150*zs+EOS050)*zt & 680 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 681 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 682 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 683 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 684 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 685 ! 686 zn_sto(jsmp) = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 687 END DO 688 ! 689 ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 690 prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 691 DO jsmp = 1, nn_sto_eos*2 692 prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp) ! potential density referenced at the surface 693 ! 694 prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_rho0 - 1._wp ) ! density anomaly (masked) 695 END DO 696 prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 697 prd (ji,jj,jk) = 0.5_wp * prd (ji,jj,jk) * ztm / nn_sto_eos 698 END_3D 699 DEALLOCATE(zn0_sto,zn_sto,zsign) 700 ! Non-stochastic equation of state 701 ELSE 702 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 703 ! 704 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 705 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 706 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 707 ztm = tmask(ji,jj,jk) ! tmask 708 ! 709 zn3 = EOS013*zt & 710 & + EOS103*zs+EOS003 711 ! 712 zn2 = (EOS022*zt & 713 & + EOS112*zs+EOS012)*zt & 714 & + (EOS202*zs+EOS102)*zs+EOS002 715 ! 716 zn1 = (((EOS041*zt & 717 & + EOS131*zs+EOS031)*zt & 718 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 719 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 720 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 721 ! 722 zn0 = (((((EOS060*zt & 723 & + EOS150*zs+EOS050)*zt & 724 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 725 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 726 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 727 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 728 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 729 ! 730 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 731 ! 732 prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface 733 ! 734 prd(ji,jj,jk) = ( zn * r1_rho0 - 1._wp ) * ztm ! density anomaly (masked) 735 END_3D 736 ENDIF 737 738 CASE( np_seos ) !== simplified EOS ==! 739 ! 740 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 741 zt = pts (ji,jj,jk,jp_tem) - 10._wp 742 zs = pts (ji,jj,jk,jp_sal) - 35._wp 743 zh = pdep (ji,jj,jk) 744 ztm = tmask(ji,jj,jk) 745 ! ! potential density referenced at the surface 746 zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt & 747 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs & 748 & - rn_nu * zt * zs 749 prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 750 ! ! density anomaly (masked) 751 zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 752 prd(ji,jj,jk) = zn * r1_rho0 * ztm 753 ! 754 END_3D 755 ! 756 END SELECT 757 ! 758 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=REAL(prd, wp), clinfo1=' eos-pot: ', & 759 & tab3d_2=REAL(prhop, wp), clinfo2=' pot : ', kdim=jpk ) 760 ! 761 IF( ln_timing ) CALL timing_stop('eos-pot') 762 ! 763 END SUBROUTINE eos_insitu_pot_t_mixed 478 764 479 765 … … 661 947 662 948 663 SUBROUTINE rab_3d ( pts, pab, Kmm )949 SUBROUTINE rab_3d_wp( pts, pab, Kmm ) 664 950 !! 665 951 INTEGER , INTENT(in ) :: Kmm ! time level index … … 667 953 REAL(wp), DIMENSION(:,:,:,:), INTENT( out) :: pab ! thermal/haline expansion ratio 668 954 !! 669 CALL rab_3d_t( pts, is_tile(pts), pab, is_tile(pab), Kmm ) 670 END SUBROUTINE rab_3d 671 672 673 SUBROUTINE rab_3d_t( pts, ktts, pab, ktab, Kmm ) 955 CALL rab_3d_t_wp( pts, is_tile(pts), pab, is_tile(pab), Kmm ) 956 END SUBROUTINE rab_3d_wp 957 958 SUBROUTINE rab_3d_mixed( pts, pab, Kmm ) 959 !! 960 INTEGER , INTENT(in ) :: Kmm ! time level index 961 REAL(dp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! pot. temperature & salinity 962 REAL(sp), DIMENSION(:,:,:,:), INTENT( out) :: pab ! thermal/haline expansion ratio 963 !! 964 CALL rab_3d_t_mixed( pts, is_tile(pts), pab, is_tile(pab), Kmm ) 965 END SUBROUTINE rab_3d_mixed 966 967 968 SUBROUTINE rab_3d_t_wp( pts, ktts, pab, ktab, Kmm ) 674 969 !!---------------------------------------------------------------------- 675 970 !! *** ROUTINE rab_3d *** … … 775 1070 IF( ln_timing ) CALL timing_stop('rab_3d') 776 1071 ! 777 END SUBROUTINE rab_3d_t 778 779 780 SUBROUTINE rab_2d( pts, pdep, pab, Kmm ) 781 !! 782 INTEGER , INTENT(in ) :: Kmm ! time level index 783 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pts ! pot. temperature & salinity 784 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pdep ! depth [m] 785 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pab ! thermal/haline expansion ratio 786 !! 787 CALL rab_2d_t(pts, is_tile(pts), pdep, is_tile(pdep), pab, is_tile(pab), Kmm) 788 END SUBROUTINE rab_2d 789 790 791 SUBROUTINE rab_2d_t( pts, ktts, pdep, ktdep, pab, ktab, Kmm ) 792 !!---------------------------------------------------------------------- 793 !! *** ROUTINE rab_2d *** 794 !! 795 !! ** Purpose : Calculates thermal/haline expansion ratio for a 2d field (unmasked) 1072 END SUBROUTINE rab_3d_t_wp 1073 1074 SUBROUTINE rab_3d_t_mixed( pts, ktts, pab, ktab, Kmm ) 1075 !!---------------------------------------------------------------------- 1076 !! *** ROUTINE rab_3d *** 1077 !! 1078 !! ** Purpose : Calculates thermal/haline expansion ratio at T-points 1079 !! 1080 !! ** Method : calculates alpha / beta at T-points 796 1081 !! 797 1082 !! ** Action : - pab : thermal/haline expansion ratio at T-points 798 1083 !!---------------------------------------------------------------------- 799 INTEGER , INTENT(in ) :: Kmm ! time level index 800 INTEGER , INTENT(in ) :: ktts, ktdep, ktab 801 REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in ) :: pts ! pot. temperature & salinity 802 REAL(wp), DIMENSION(A2D_T(ktdep) ), INTENT(in ) :: pdep ! depth [m] 803 REAL(wp), DIMENSION(A2D_T(ktab),JPTS), INTENT( out) :: pab ! thermal/haline expansion ratio 1084 INTEGER , INTENT(in ) :: Kmm ! time level index 1085 INTEGER , INTENT(in ) :: ktts, ktab 1086 REAL(dp), DIMENSION(A2D_T(ktts),JPK,JPTS), INTENT(in ) :: pts ! pot. temperature & salinity 1087 REAL(sp), DIMENSION(A2D_T(ktab),JPK,JPTS), INTENT( out) :: pab ! thermal/haline expansion ratio 804 1088 ! 805 1089 INTEGER :: ji, jj, jk ! dummy loop indices 806 REAL(wp) :: zt , zh , zs 1090 REAL(wp) :: zt , zh , zs , ztm ! local scalars 807 1091 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 808 1092 !!---------------------------------------------------------------------- 809 1093 ! 810 IF( ln_timing ) CALL timing_start('rab_2d') 811 ! 812 pab(:,:,:) = 0._wp 1094 IF( ln_timing ) CALL timing_start('rab_3d') 813 1095 ! 814 1096 SELECT CASE ( neos ) … … 816 1098 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 817 1099 ! 818 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 819 ! 820 zh = pdep(ji,jj) * r1_Z0 ! depth 821 zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature 822 zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 1100 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 1101 ! 1102 zh = gdept(ji,jj,jk,Kmm) * r1_Z0 ! depth 1103 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 1104 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 1105 ztm = tmask(ji,jj,jk) ! tmask 823 1106 ! 824 1107 ! alpha … … 841 1124 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 842 1125 ! 843 pab(ji,jj,j p_tem) = zn * r1_rho01126 pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm 844 1127 ! 845 1128 ! beta … … 862 1145 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 863 1146 ! 1147 pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm 1148 ! 1149 END_3D 1150 ! 1151 CASE( np_seos ) !== simplified EOS ==! 1152 ! 1153 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 1154 zt = pts (ji,jj,jk,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) 1155 zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) 1156 zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point 1157 ztm = tmask(ji,jj,jk) ! land/sea bottom mask = surf. mask 1158 ! 1159 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 1160 pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm ! alpha 1161 ! 1162 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 1163 pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm ! beta 1164 ! 1165 END_3D 1166 ! 1167 CASE DEFAULT 1168 WRITE(ctmp1,*) ' bad flag value for neos = ', neos 1169 CALL ctl_stop( 'rab_3d:', ctmp1 ) 1170 ! 1171 END SELECT 1172 ! 1173 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=REAL(pab(:,:,:,jp_tem), wp), clinfo1=' rab_3d_t: ', & 1174 & tab3d_2=REAL(pab(:,:,:,jp_sal), wp), clinfo2=' rab_3d_s : ', kdim=jpk ) 1175 ! 1176 IF( ln_timing ) CALL timing_stop('rab_3d') 1177 ! 1178 END SUBROUTINE rab_3d_t_mixed 1179 1180 1181 SUBROUTINE rab_2d( pts, pdep, pab, Kmm ) 1182 !! 1183 INTEGER , INTENT(in ) :: Kmm ! time level index 1184 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pts ! pot. temperature & salinity 1185 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pdep ! depth [m] 1186 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pab ! thermal/haline expansion ratio 1187 !! 1188 CALL rab_2d_t(pts, is_tile(pts), pdep, is_tile(pdep), pab, is_tile(pab), Kmm) 1189 END SUBROUTINE rab_2d 1190 1191 1192 SUBROUTINE rab_2d_t( pts, ktts, pdep, ktdep, pab, ktab, Kmm ) 1193 !!---------------------------------------------------------------------- 1194 !! *** ROUTINE rab_2d *** 1195 !! 1196 !! ** Purpose : Calculates thermal/haline expansion ratio for a 2d field (unmasked) 1197 !! 1198 !! ** Action : - pab : thermal/haline expansion ratio at T-points 1199 !!---------------------------------------------------------------------- 1200 INTEGER , INTENT(in ) :: Kmm ! time level index 1201 INTEGER , INTENT(in ) :: ktts, ktdep, ktab 1202 REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in ) :: pts ! pot. temperature & salinity 1203 REAL(wp), DIMENSION(A2D_T(ktdep) ), INTENT(in ) :: pdep ! depth [m] 1204 REAL(wp), DIMENSION(A2D_T(ktab),JPTS), INTENT( out) :: pab ! thermal/haline expansion ratio 1205 ! 1206 INTEGER :: ji, jj, jk ! dummy loop indices 1207 REAL(wp) :: zt , zh , zs ! local scalars 1208 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 1209 !!---------------------------------------------------------------------- 1210 ! 1211 IF( ln_timing ) CALL timing_start('rab_2d') 1212 ! 1213 pab(:,:,:) = 0._wp 1214 ! 1215 SELECT CASE ( neos ) 1216 ! 1217 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 1218 ! 1219 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1220 ! 1221 zh = pdep(ji,jj) * r1_Z0 ! depth 1222 zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature 1223 zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 1224 ! 1225 ! alpha 1226 zn3 = ALP003 1227 ! 1228 zn2 = ALP012*zt + ALP102*zs+ALP002 1229 ! 1230 zn1 = ((ALP031*zt & 1231 & + ALP121*zs+ALP021)*zt & 1232 & + (ALP211*zs+ALP111)*zs+ALP011)*zt & 1233 & + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 1234 ! 1235 zn0 = ((((ALP050*zt & 1236 & + ALP140*zs+ALP040)*zt & 1237 & + (ALP230*zs+ALP130)*zs+ALP030)*zt & 1238 & + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt & 1239 & + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt & 1240 & + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 1241 ! 1242 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 1243 ! 1244 pab(ji,jj,jp_tem) = zn * r1_rho0 1245 ! 1246 ! beta 1247 zn3 = BET003 1248 ! 1249 zn2 = BET012*zt + BET102*zs+BET002 1250 ! 1251 zn1 = ((BET031*zt & 1252 & + BET121*zs+BET021)*zt & 1253 & + (BET211*zs+BET111)*zs+BET011)*zt & 1254 & + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 1255 ! 1256 zn0 = ((((BET050*zt & 1257 & + BET140*zs+BET040)*zt & 1258 & + (BET230*zs+BET130)*zs+BET030)*zt & 1259 & + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt & 1260 & + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt & 1261 & + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 1262 ! 1263 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 1264 ! 864 1265 pab(ji,jj,jp_sal) = zn / zs * r1_rho0 865 1266 ! … … 997 1398 !! 998 1399 INTEGER , INTENT(in ) :: Kmm ! time level index 999 REAL( wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature and salinity [Celsius,psu]1400 REAL(dp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature and salinity [Celsius,psu] 1000 1401 REAL(wp), DIMENSION(:,:,:,:) , INTENT(in ) :: pab ! thermal/haline expansion coef. [Celsius-1,psu-1] 1001 1402 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: pn2 ! Brunt-Vaisala frequency squared [1/s^2] … … 1021 1422 INTEGER , INTENT(in ) :: Kmm ! time level index 1022 1423 INTEGER , INTENT(in ) :: ktab, ktn2 1023 REAL( wp), DIMENSION(jpi,jpj, jpk,jpts), INTENT(in ) :: pts ! pot. temperature and salinity [Celsius,psu]1424 REAL(dp), DIMENSION(jpi,jpj, jpk,jpts), INTENT(in ) :: pts ! pot. temperature and salinity [Celsius,psu] 1024 1425 REAL(wp), DIMENSION(A2D_T(ktab),JPK,JPTS), INTENT(in ) :: pab ! thermal/haline expansion coef. [Celsius-1,psu-1] 1025 1426 REAL(wp), DIMENSION(A2D_T(ktn2),JPK ), INTENT( out) :: pn2 ! Brunt-Vaisala frequency squared [1/s^2] -
NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/tranpc.F90
r14986 r15066 101 101 ! 102 102 CALL eos_rab( CASTWP(pts(:,:,:,:,Kaa)), zab, Kmm ) ! after alpha and beta (given on T-points) 103 CALL bn2 ( CASTWP(pts(:,:,:,:,Kaa)), zab, zn2, Kmm ) ! after Brunt-Vaisala (given on W-points)103 CALL bn2 ( pts(:,:,:,:,Kaa), zab, zn2, Kmm ) ! after Brunt-Vaisala (given on W-points) 104 104 ! 105 105 IF( .NOT. l_istiled .OR. ntile == 1 ) nnpcc = 0 ! Do only on the first tile -
NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/ZDF/zdfosm.F90
r14986 r15066 3295 3295 ! w-level of the mixing and mixed layers 3296 3296 CALL eos_rab( CASTWP(ts(:,:,:,:,Kmm)), rab_n, Kmm ) 3297 CALL bn2( CASTWP(ts(:,:,:,:,Kmm)), rab_n, rn2, Kmm )3297 CALL bn2( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm ) 3298 3298 imld_rst(:,:) = nlb10 ! Initialization to the number of w ocean point 3299 3299 hbl(:,:) = 0.0_wp ! Here hbl used as a dummy variable, integrating vertically N^2 -
NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/step.F90
r14986 r15066 171 171 CALL eos_rab( CASTWP(ts(:,:,:,:,Nbb)), rab_b, Nnn ) ! before local thermal/haline expension ratio at T-points 172 172 CALL eos_rab( CASTWP(ts(:,:,:,:,Nnn)), rab_n, Nnn ) ! now local thermal/haline expension ratio at T-points 173 CALL bn2 ( CASTWP(ts(:,:,:,:,Nbb)), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency174 CALL bn2 ( CASTWP(ts(:,:,:,:,Nnn)), rab_n, rn2, Nnn ) ! now Brunt-Vaisala frequency173 CALL bn2 ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency 174 CALL bn2 ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn ) ! now Brunt-Vaisala frequency 175 175 176 176 ! VERTICAL PHYSICS -
NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/stpmlf.F90
r14986 r15066 176 176 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 177 177 ! THERMODYNAMICS 178 CALL eos_rab( CASTWP(ts(:,:,:,:,Nbb)), rab_b, Nnn ) ! before local thermal/haline expension ratio at T-points179 CALL eos_rab( CASTWP(ts(:,:,:,:,Nnn)), rab_n, Nnn ) ! now local thermal/haline expension ratio at T-points180 CALL bn2 ( CASTWP(ts(:,:,:,:,Nbb)), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency181 CALL bn2 ( CASTWP(ts(:,:,:,:,Nnn)), rab_n, rn2, Nnn ) ! now Brunt-Vaisala frequency178 CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn ) ! before local thermal/haline expension ratio at T-points 179 CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn ) ! now local thermal/haline expension ratio at T-points 180 CALL bn2 ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency 181 CALL bn2 ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn ) ! now Brunt-Vaisala frequency 182 182 183 183 ! VERTICAL PHYSICS … … 217 217 CALL ssh_nxt ( kstp, Nbb, Nnn, ssh, Naa ) ! after ssh (includes call to div_hor) 218 218 IF( .NOT.lk_linssh ) THEN 219 CALL dom_qco_r3c( CASTWP(ssh(:,:,Naa)), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa) ) ! "after" ssh/h_0 ratio at t,u,v pts219 CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa) ) ! "after" ssh/h_0 ratio at t,u,v pts 220 220 IF( ln_dynspg_exp ) & 221 & CALL dom_qco_r3c( CASTWP(ssh(:,:,Nnn)), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn), r3f(:,:) ) ! spg_exp : needed only for "now" ssh/h_0 ratio at f point221 & CALL dom_qco_r3c( ssh(:,:,Nnn), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn), r3f(:,:) ) ! spg_exp : needed only for "now" ssh/h_0 ratio at f point 222 222 ENDIF 223 223 CALL wzv ( kstp, Nbb, Nnn, Naa, ww ) ! Nnn cross-level velocity … … 227 227 zgdept(:,:,jk) = gdept(:,:,jk,Nnn) 228 228 END DO 229 CALL eos ( CASTWP(ts(:,:,:,:,Nnn)), rhd, rhop, zgdept ) ! now in situ density for hpg computation229 CALL eos ( ts(:,:,:,:,Nnn), rhd, rhop, zgdept ) ! now in situ density for hpg computation 230 230 DEALLOCATE( zgdept ) 231 231 … … 268 268 ! as well as vertical scale factors and vertical velocity need to be updated 269 269 CALL div_hor ( kstp, Nbb, Nnn ) ! Horizontal divergence (2nd call in time-split case) 270 IF(.NOT.lk_linssh) CALL dom_qco_r3c( CASTWP(ssh(:,:,Naa)), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) ) ! update ssh/h_0 ratio at t,u,v,f pts270 IF(.NOT.lk_linssh) CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) ) ! update ssh/h_0 ratio at t,u,v,f pts 271 271 ENDIF 272 272 CALL dyn_zdf ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa ) ! vertical diffusion … … 303 303 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 304 304 CALL ssh_atf ( kstp, Nbb, Nnn, Naa, ssh ) ! time filtering of "now" sea surface height 305 IF(.NOT.lk_linssh) CALL dom_qco_r3c( CASTWP(ssh(:,:,Nnn)), r3t_f, r3u_f, r3v_f ) ! "now" ssh/h_0 ratio from filtrered ssh305 IF(.NOT.lk_linssh) CALL dom_qco_r3c( ssh(:,:,Nnn), r3t_f, r3u_f, r3v_f ) ! "now" ssh/h_0 ratio from filtrered ssh 306 306 #if defined key_top 307 307 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Note: See TracChangeset
for help on using the changeset viewer.
