Changeset 3618 for branches/2012/dev_UKMO_2012/NEMOGCM/NEMO
- Timestamp:
- 2012-11-20T19:20:57+01:00 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r3600 r3618 15 15 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 16 16 !! 3.3 ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level 17 !! 3.4 ! 2012-08 (J. Siddorn) added Siddorn a dn Furner stretching functio17 !! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function 18 18 !!---------------------------------------------------------------------- 19 19 … … 28 28 !! zgr_zps : z-coordinate with partial steps 29 29 !! zgr_sco : s-coordinate 30 !! fssig : sigma coordinate non-dimensionalfunction31 !! f gamma : Siddorn and Furner stretchingfunction32 !! dfssig : derivative of the sigma coordinate function !!gm (currently missing!)30 !! fssig : tanh stretch function 31 !! fssig1 : Song and Haidvogel 1994 stretch function 32 !! fgamma : Siddorn and Furner 2012 stretching function 33 33 !!--------------------------------------------------------------------- 34 34 USE oce ! ocean variables … … 51 51 LOGICAL :: ln_s_sh94 = .false. ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T) 52 52 LOGICAL :: ln_s_sf12 = .true. ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T) 53 LOGICAL :: ln_sigcrit = .false. ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch54 53 ! 55 54 REAL(wp) :: rn_sbot_min = 300._wp ! minimum depth of s-bottom surface (>0) (m) … … 63 62 ! ! ( rn_bb=0; top only, rn_bb =1; top and bottom) 64 63 ! Siddorn and Furner stretching parameters 64 LOGICAL :: ln_sigcrit = .false. ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 65 65 REAL(wp) :: rn_alpha = 4.4_wp ! control parameter ( > 1 stretch towards surface, < 1 towards seabed) 66 66 REAL(wp) :: rn_efold = 0.0_wp ! efold length scale for transition to stretched coord … … 1065 1065 !! hbatv = mj( hbatt ) 1066 1066 !! hbatf = mi( mj( hbatt ) ) 1067 !! - Compute gsigt, gsigw, esigt,esigw from an analytical1067 !! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical 1068 1068 !! function and its derivative given as function. 1069 !! gsigt(k) = fssig (k )1070 !! gsigw(k) = fssig (k-0.5)1071 !! esigt(k) = fsdsig(k )1072 !! esigw(k) = fsdsig(k-0.5)1069 !! z_gsigt(k) = fssig (k ) 1070 !! z_gsigw(k) = fssig (k-0.5) 1071 !! z_esigt(k) = fsdsig(k ) 1072 !! z_esigw(k) = fsdsig(k-0.5) 1073 1073 !! Three options for stretching are give, and they can be modified 1074 1074 !! following the users requirements. Nevertheless, the output as … … 1125 1125 WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb 1126 1126 WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12 1127 WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit 1127 1128 WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients' 1128 1129 WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha … … 1383 1384 & ' MAX ', MAXVAL( mbathy(:,:) ) 1384 1385 1385 ! ! =============1386 IF(lwp) THEN ! Control print1387 ! ! =============1388 WRITE(numout,*)1389 WRITE(numout,*) ' domzgr: vertical coefficients for model level'1390 WRITE(numout, "(9x,' level gsigt gsigw esigt esigw gsi3w')" )1391 WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk )1392 ENDIF1393 1386 IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain 1394 1387 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) … … 1508 1501 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 1509 1502 ! 1510 REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3,gsi3w31511 REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3,esigwv31512 1513 CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3,gsi3w3 )1514 CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3,esigwv3 )1515 1516 gsigw3 = 0._wp ; gsigt3 = 0._wp ;gsi3w3 = 0._wp1517 esigt3 = 0._wp ;esigw3 = 0._wp1518 esigtu3 = 0._wp ; esigtv3 = 0._wp ;esigtf3 = 0._wp1519 esigwu3 = 0._wp ;esigwv3 = 0._wp1503 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 1504 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 1505 1506 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 1507 CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 1508 1509 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp 1510 z_esigt3 = 0._wp ; z_esigw3 = 0._wp 1511 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp 1512 z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp 1520 1513 1521 1514 DO ji = 1, jpi … … 1524 1517 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 1525 1518 DO jk = 1, jpk 1526 gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )1527 gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )1519 z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 1520 z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) 1528 1521 END DO 1529 1522 ELSE ! shallow water, uniform sigma 1530 1523 DO jk = 1, jpk 1531 gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp)1532 gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)1524 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp) 1525 z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 1533 1526 END DO 1534 1527 ENDIF 1535 1528 ! 1536 1529 DO jk = 1, jpkm1 1537 esigt3(ji,jj,jk ) = gsigw3(ji,jj,jk+1) -gsigw3(ji,jj,jk)1538 esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) -gsigt3(ji,jj,jk)1539 END DO 1540 esigw3(ji,jj,1 ) = 2._wp * ( gsigt3(ji,jj,1 ) -gsigw3(ji,jj,1 ) )1541 esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) -gsigw3(ji,jj,jpk) )1530 z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk) 1531 z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 1532 END DO 1533 z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) ) 1534 z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) ) 1542 1535 ! 1543 1536 ! Coefficients for vertical depth as the sum of e3w scale factors 1544 gsi3w3(ji,jj,1) = 0.5_wp *esigw3(ji,jj,1)1537 z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1) 1545 1538 DO jk = 2, jpk 1546 gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) +esigw3(ji,jj,jk)1539 z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 1547 1540 END DO 1548 1541 ! … … 1550 1543 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1551 1544 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1552 gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)* gsigt3(ji,jj,jk)+rn_hc*zcoeft )1553 gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)* gsigw3(ji,jj,jk)+rn_hc*zcoefw )1554 gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)* gsi3w3(ji,jj,jk)+rn_hc*zcoeft )1545 gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 1546 gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 1547 gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 1555 1548 END DO 1556 1549 ! … … 1561 1554 DO jj = 1, jpjm1 1562 1555 DO jk = 1, jpk 1563 esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) &1556 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 1564 1557 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1565 esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) &1558 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 1566 1559 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1567 esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) &1568 & + hbatt(ji,jj+1)* esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) &1560 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & 1561 & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 1569 1562 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1570 esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) &1563 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 1571 1564 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1572 esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) &1565 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 1573 1566 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1574 1567 ! 1575 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)* esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )1576 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)* esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )1577 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)* esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )1578 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)* esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )1568 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1569 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1570 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1571 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1579 1572 ! 1580 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)* esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )1581 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)* esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )1582 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)* esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )1573 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1574 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1575 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1583 1576 END DO 1584 1577 END DO 1585 1578 END DO 1586 1579 1587 CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3,gsi3w3 )1588 CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3,esigwv3 )1580 CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 1581 CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 1589 1582 1590 1583 END SUBROUTINE s_sh94 … … 1609 1602 ! 1610 1603 INTEGER :: ji, jj, jk ! dummy loop argument 1611 REAL(wp) :: fsmth! smoothing around critical depth1612 REAL(wp) :: z ss, zbb! Surface and bottom cell thickness in sigma space1613 ! 1614 REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3,gsi3w31615 REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3,esigwv31616 1617 ! 1618 CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3,gsi3w3 )1619 CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3,esigwv3 )1620 1621 gsigw3 = 0._wp ; gsigt3 = 0._wp ;gsi3w3 = 0._wp1622 esigt3 = 0._wp ;esigw3 = 0._wp1623 esigtu3 = 0._wp ; esigtv3 = 0._wp ;esigtf3 = 0._wp1624 esigwu3 = 0._wp ;esigwv3 = 0._wp1604 REAL(wp) :: zsmth ! smoothing around critical depth 1605 REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space 1606 ! 1607 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 1608 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 1609 1610 ! 1611 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 1612 CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 1613 1614 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp 1615 z_esigt3 = 0._wp ; z_esigw3 = 0._wp 1616 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp 1617 z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp 1625 1618 1626 1619 DO ji = 1, jpi … … 1629 1622 IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma 1630 1623 1631 z bb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,.1624 zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,. 1632 1625 ! could be changed by users but care must be taken to do so carefully 1633 z bb = 1.0_wp-(zbb/hbatt(ji,jj))1626 zzb = 1.0_wp-(zzb/hbatt(ji,jj)) 1634 1627 1635 z ss = rn_zs / hbatt(ji,jj)1628 zzs = rn_zs / hbatt(ji,jj) 1636 1629 1637 1630 IF (rn_efold /= 0.0_wp) THEN 1638 fsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )1631 zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) 1639 1632 ELSE 1640 fsmth = 1.0_wp1633 zsmth = 1.0_wp 1641 1634 ENDIF 1642 1635 1643 1636 DO jk = 1, jpk 1644 gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)1645 gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)1637 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) 1638 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp) 1646 1639 ENDDO 1647 gsigw3(ji,jj,:) = fgamma( gsigw3(ji,jj,:), zbb, zss, fsmth )1648 gsigt3(ji,jj,:) = fgamma( gsigt3(ji,jj,:), zbb, zss, fsmth )1640 z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth ) 1641 z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth ) 1649 1642 1650 1643 ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma 1651 1644 1652 1645 DO jk = 1, jpk 1653 gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)1654 gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)1646 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) 1647 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 1655 1648 END DO 1656 1649 … … 1658 1651 1659 1652 DO jk = 1, jpk 1660 gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))1661 gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))1653 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 1654 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 1662 1655 END DO 1663 1656 … … 1665 1658 1666 1659 DO jk = 1, jpkm1 1667 esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) -gsigw3(ji,jj,jk)1668 esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) -gsigt3(ji,jj,jk)1660 z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk) 1661 z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 1669 1662 END DO 1670 esigw3(ji,jj,1 ) = 2.0_wp * (gsigt3(ji,jj,1 ) -gsigw3(ji,jj,1 ))1671 esigt3(ji,jj,jpk) = 2.0_wp * (gsigt3(ji,jj,jpk) -gsigw3(ji,jj,jpk))1663 z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 )) 1664 z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk)) 1672 1665 1673 1666 ! Coefficients for vertical depth as the sum of e3w scale factors 1674 gsi3w3(ji,jj,1) = 0.5 *esigw3(ji,jj,1)1667 z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1) 1675 1668 DO jk = 2, jpk 1676 gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) +esigw3(ji,jj,jk)1669 z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 1677 1670 END DO 1678 1671 1679 1672 DO jk = 1, jpk 1680 gdept (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))* gsigt3(ji,jj,jk)1681 gdepw (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))* gsigw3(ji,jj,jk)1682 gdep3w(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))* gsi3w3(ji,jj,jk)1673 gdept (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 1674 gdepw (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 1675 gdep3w(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 1683 1676 END DO 1684 1677 … … 1690 1683 1691 1684 DO jk = 1, jpk 1692 esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) / &1685 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 1693 1686 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1694 esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) / &1687 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 1695 1688 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1696 esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) + &1697 hbatt(ji,jj+1)* esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) / &1689 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) + & 1690 hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 1698 1691 ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1699 esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) / &1692 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 1700 1693 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1701 esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) / &1694 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 1702 1695 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1703 1696 1704 e3t(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))* esigt3(ji,jj,jk)1705 e3u(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))* esigtu3(ji,jj,jk)1706 e3v(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))* esigtv3(ji,jj,jk)1707 e3f(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))* esigtf3(ji,jj,jk)1697 e3t(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) 1698 e3u(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk) 1699 e3v(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk) 1700 e3f(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk) 1708 1701 ! 1709 e3w(ji,jj,jk)=hbatt(ji,jj)* esigw3(ji,jj,jk)1710 e3uw(ji,jj,jk)=hbatu(ji,jj)* esigwu3(ji,jj,jk)1711 e3vw(ji,jj,jk)=hbatv(ji,jj)* esigwv3(ji,jj,jk)1702 e3w(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk) 1703 e3uw(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk) 1704 e3vw(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk) 1712 1705 END DO 1713 1706 1714 1707 ENDDO 1715 1708 ENDDO 1716 1717 CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3 ) 1718 CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 1709 ! ! ============= 1710 1711 CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 1712 CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 1719 1713 1720 1714 END SUBROUTINE s_sf12 … … 1735 1729 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 1736 1730 1731 REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w 1732 REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw 1733 1734 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 1735 CALL wrk_alloc( jpk, z_esigt, z_esigw ) 1736 1737 z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp 1738 z_esigt = 0._wp ; z_esigw = 0._wp 1739 1737 1740 DO jk = 1, jpk 1738 gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )1739 gsigt(jk) = -fssig( REAL(jk,wp) )1740 END DO 1741 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' gsigw 1 jpk ', gsigw(1),gsigw(jpk)1741 z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 1742 z_gsigt(jk) = -fssig( REAL(jk,wp) ) 1743 END DO 1744 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk) 1742 1745 ! 1743 1746 ! Coefficients for vertical scale factors at w-, t- levels … … 1745 1748 !!gm or betteroffer the 2 possibilities.... 1746 1749 DO jk = 1, jpkm1 1747 esigt(jk ) = gsigw(jk+1) -gsigw(jk)1748 esigw(jk+1) = gsigt(jk+1) -gsigt(jk)1749 END DO 1750 esigw( 1 ) = 2._wp * ( gsigt(1 ) -gsigw(1 ) )1751 esigt(jpk) = 2._wp * ( gsigt(jpk) -gsigw(jpk) )1750 z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk) 1751 z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk) 1752 END DO 1753 z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) ) 1754 z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) ) 1752 1755 ! 1753 1756 ! Coefficients for vertical depth as the sum of e3w scale factors 1754 gsi3w(1) = 0.5_wp *esigw(1)1757 z_gsi3w(1) = 0.5_wp * z_esigw(1) 1755 1758 DO jk = 2, jpk 1756 gsi3w(jk) = gsi3w(jk-1) +esigw(jk)1759 z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk) 1757 1760 END DO 1758 1761 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) … … 1760 1763 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1761 1764 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1762 gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))* gsigt(jk) + hift(:,:)*zcoeft )1763 gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))* gsigw(jk) + hift(:,:)*zcoefw )1764 gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))* gsi3w(jk) + hift(:,:)*zcoeft )1765 gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 1766 gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 1767 gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 1765 1768 END DO 1766 1769 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) … … 1768 1771 DO ji = 1, jpi 1769 1772 DO jk = 1, jpk 1770 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))* esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )1771 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))* esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )1772 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))* esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )1773 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))* esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )1773 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1774 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1775 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1776 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 1774 1777 ! 1775 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1776 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1777 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1778 END DO 1779 END DO 1780 END DO 1778 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1779 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1780 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1781 END DO 1782 END DO 1783 END DO 1784 1785 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 1786 CALL wrk_dealloc( jpk, z_esigt, z_esigw ) 1787 1781 1788 END SUBROUTINE s_tanh 1782 1789 … … 1832 1839 1833 1840 1834 FUNCTION fgamma( pk1, Zbb, Zss, F ) RESULT( gam)1841 FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma ) 1835 1842 !!---------------------------------------------------------------------- 1836 1843 !! *** ROUTINE fgamma *** … … 1849 1856 !! Reference : Siddorn and Furner, in prep 1850 1857 !!---------------------------------------------------------------------- 1851 REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate1852 REAL(wp) :: gam(jpk) ! stretched coordinate1853 REAL(wp), INTENT(in ) :: Zbb! Bottom box depth1854 REAL(wp), INTENT(in ) :: Zss! surface box depth1855 REAL(wp), INTENT(in ) :: F! Smoothing parameter1856 REAL(wp) :: a1,a2,a3! local variables1857 REAL(wp) :: n1,n2! local variables1858 REAL(wp) :: A,B,X! local variables1858 REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate 1859 REAL(wp) :: p_gamma(jpk) ! stretched coordinate 1860 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth 1861 REAL(wp), INTENT(in ) :: pzs ! surface box depth 1862 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter 1863 REAL(wp) :: za1,za2,za3 ! local variables 1864 REAL(wp) :: zn1,zn2 ! local variables 1865 REAL(wp) :: za,zb,zx ! local variables 1859 1866 integer :: jk 1860 1867 !!---------------------------------------------------------------------- 1861 1868 ! 1862 1869 1863 n1 = 1./(jpk-1.)1864 n2 = 1. -n11865 1866 a1 = (rn_alpha+2.0_wp)*n1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n1**(rn_alpha+2.0_wp)1867 a2 = (rn_alpha+2.0_wp)*n2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n2**(rn_alpha+2.0_wp)1868 a3 = ( n2**3.0_wp - a2)/( n1**3.0_wp -a1)1870 zn1 = 1./(jpk-1.) 1871 zn2 = 1. - zn1 1872 1873 za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 1874 za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp) 1875 za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) 1869 1876 1870 A = Zbb - a3*(Zss-a1)-a21871 A = A/( n2-0.5_wp*(a2+n2**2.0_wp) - a3*(n1-0.5_wp*(a1+n1**2.0_wp) ) )1872 B = (Zss - a1 - A*( n1-0.5_wp*(a1+n1**2.0_wp ) ) ) / (n1**3.0_wp -a1)1873 X = 1.0_wp-A/2.0_wp-B1877 za = pzb - za3*(pzs-za1)-za2 1878 za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) ) 1879 zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) 1880 zx = 1.0_wp-za/2.0_wp-zb 1874 1881 1875 1882 DO jk = 1, jpk 1876 gam(jk) = A*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+B*pk1(jk)**3.0_wp + X*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )1877 gam(jk) = gam(jk)*F+pk1(jk)*(1.0_wp-F)1883 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) ) 1884 p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth) 1878 1885 ENDDO 1879 1886
Note: See TracChangeset
for help on using the changeset viewer.
