[3] | 1 | MODULE diahdy |
---|
| 2 | !!====================================================================== |
---|
| 3 | !! *** MODULE diahdy *** |
---|
| 4 | !! Ocean diagnostics : computation the dynamical heigh |
---|
| 5 | !!====================================================================== |
---|
| 6 | #if defined key_diahdy || defined key_esopa |
---|
| 7 | !!---------------------------------------------------------------------- |
---|
| 8 | !! 'key_diahdy' : dynamical heigh diagnostics |
---|
| 9 | !!---------------------------------------------------------------------- |
---|
| 10 | !! dia_hdy : dynamical heigh computation |
---|
| 11 | !!---------------------------------------------------------------------- |
---|
| 12 | !! * Modules used |
---|
| 13 | USE oce ! ocean dynamics and tracers |
---|
| 14 | USE dom_oce ! ocean space and time domain |
---|
| 15 | USE phycst ! physical constants |
---|
| 16 | USE in_out_manager ! I/O manager |
---|
| 17 | |
---|
| 18 | IMPLICIT NONE |
---|
[32] | 19 | PRIVATE |
---|
[3] | 20 | |
---|
[32] | 21 | !! * Routine accessibility |
---|
| 22 | PUBLIC dia_hdy ! called in step.F90 module |
---|
| 23 | |
---|
[3] | 24 | !! * Shared module variables |
---|
[32] | 25 | LOGICAL, PUBLIC, PARAMETER :: lk_diahdy = .TRUE. !: dynamical heigh flag |
---|
[3] | 26 | |
---|
| 27 | !! * Module variables |
---|
| 28 | REAL(wp), DIMENSION(jpk) :: & |
---|
| 29 | rhosp ! ??? |
---|
| 30 | |
---|
| 31 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: & |
---|
| 32 | hdy ! dynamical heigh |
---|
| 33 | |
---|
| 34 | !! * Substitutions |
---|
| 35 | # include "domzgr_substitute.h90" |
---|
| 36 | !!---------------------------------------------------------------------- |
---|
[247] | 37 | !! OPA 9.0 , LOCEAN-IPSL (2005) |
---|
[719] | 38 | !! $Header$ |
---|
[247] | 39 | !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt |
---|
[3] | 40 | !!---------------------------------------------------------------------- |
---|
| 41 | |
---|
| 42 | CONTAINS |
---|
| 43 | |
---|
| 44 | SUBROUTINE dia_hdy ( kt ) |
---|
| 45 | !!--------------------------------------------------------------------- |
---|
| 46 | !! *** ROUTINE dia_hdy *** |
---|
| 47 | !! |
---|
| 48 | !! ** Purpose : Computes the dynamical heigh |
---|
| 49 | !! |
---|
| 50 | !! ** Method : Millero + Poisson |
---|
| 51 | !! |
---|
| 52 | !! References : |
---|
| 53 | !! A. E. Gill, atmosphere-ocean dynamics 7.7 pp 215 |
---|
| 54 | !! |
---|
| 55 | !! History : |
---|
| 56 | !! ! 9x-xx (P. Delecluse, C. Perigaud) Original code |
---|
| 57 | !! ! 93-10 (C. Perigaud) a trapezoidal vertical integration |
---|
| 58 | !! consistent WITH the code |
---|
| 59 | !! ! 93-12 (G. Madec M. Imbard) |
---|
| 60 | !! ! 96-03 (N. Ferry) integration at t-points |
---|
| 61 | !! 8.5 ! 02-06 (G. Madec) F90: Free form and module |
---|
| 62 | !!---------------------------------------------------------------------- |
---|
| 63 | !! * Arguments |
---|
| 64 | INTEGER, INTENT( in ) :: kt ! ocean time-step index |
---|
| 65 | |
---|
| 66 | !! * Local declarations |
---|
| 67 | INTEGER :: ji, jj, jk |
---|
[32] | 68 | INTEGER :: ihdsup, ik |
---|
[3] | 69 | |
---|
| 70 | REAL(wp) :: zgdsup, za, zb, zciint, zfacto, zhd |
---|
| 71 | REAL(wp) :: zp, zh, zt, zs, zxk, zq, zsr, zr1, zr2, zr3, zr4 |
---|
[32] | 72 | REAL(wp) :: ze, zbw, zc, zd, zaw, zb1, za1, zkw, zk0 |
---|
[3] | 73 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsva |
---|
| 74 | REAL(wp), DIMENSION(jpk) :: zwkx, zwky, zwkz |
---|
| 75 | REAL(wp) :: fsatg |
---|
| 76 | REAL(wp) :: pfps, pfpt, pfphp |
---|
| 77 | |
---|
| 78 | ! Adiabatic laspse rate fsatg, defined as the change of temperature |
---|
| 79 | ! per unit pressure for adiabatic change of pressure of an element |
---|
| 80 | ! of seawater (bryden,h.,1973,deep-sea res.,20,401-408). |
---|
| 81 | ! units: |
---|
| 82 | ! pressure pfphp decibars |
---|
| 83 | ! temperature pfpt deg celsius (ipts-68) |
---|
| 84 | ! salinity pfps (ipss-78) |
---|
| 85 | ! adiabatic fsatg deg. c/decibar |
---|
| 86 | ! checkvalue: atg=3.255976e-4 c/dbar for pfps=40 (ipss-78), |
---|
| 87 | ! pfpt=40 deg c, pfphp=10000 decibars |
---|
| 88 | |
---|
| 89 | fsatg(pfps,pfpt,pfphp) & |
---|
| 90 | = (((-2.1687e-16*pfpt+1.8676e-14)*pfpt-4.6206e-13)*pfphp & |
---|
| 91 | +((2.7759e-12*pfpt-1.1351e-10)*(pfps-35.)+((-5.4481e-14*pfpt & |
---|
| 92 | +8.733e-12)*pfpt-6.7795e-10)*pfpt+1.8741e-8))*pfphp & |
---|
| 93 | +(-4.2393e-8*pfpt+1.8932e-6)*(pfps-35.) & |
---|
| 94 | +((6.6228e-10*pfpt-6.836e-8)*pfpt+8.5258e-6)*pfpt+3.5803e-5 |
---|
| 95 | !!---------------------------------------------------------------------- |
---|
| 96 | |
---|
| 97 | ! 1. height dynamic |
---|
| 98 | ! ----------------- |
---|
| 99 | ! depth for reference |
---|
| 100 | |
---|
| 101 | zgdsup = 1500. |
---|
| 102 | |
---|
| 103 | ! below for hdyn levitus |
---|
| 104 | |
---|
| 105 | IF( kt == nit000 ) THEN |
---|
| 106 | IF(lwp) WRITE(numout,*) |
---|
| 107 | IF(lwp) WRITE(numout,*) 'dia_hdy : computation of dynamical heigh' |
---|
| 108 | IF(lwp) WRITE(numout,*) '~~~~~~~' |
---|
[474] | 109 | IF( .NOT. ln_zco ) & |
---|
| 110 | & CALL ctl_stop( ' ln_zps or ln_sco, Dynamical height diagnostics not yet implemented' ) |
---|
[3] | 111 | DO jk = 1, jpk |
---|
[460] | 112 | IF( gdepw_0(jk) > zgdsup ) GOTO 110 |
---|
[3] | 113 | END DO |
---|
| 114 | IF(lwp) WRITE(numout,*)'problem zgdsup greater than gdepw(jpk)' |
---|
| 115 | STOP 'dia_hdy' |
---|
| 116 | 110 CONTINUE |
---|
| 117 | ihdsup = jk - 1 |
---|
| 118 | IF(lwp) WRITE(numout,*)' ihdsup = ', ihdsup |
---|
| 119 | |
---|
| 120 | ! Interpolation coefficients for zgdsup-gdepw(ihdsup) layer |
---|
| 121 | |
---|
[460] | 122 | za = gdepw_0(ihdsup ) |
---|
| 123 | zb = gdepw_0(ihdsup+1) |
---|
[3] | 124 | IF( za > zgdsup .OR. zb < zgdsup ) THEN |
---|
| 125 | IF(lwp) WRITE(numout,*) za, zb, ihdsup, zgdsup |
---|
| 126 | IF(lwp) WRITE(numout,*) ' bad ihdsup' |
---|
| 127 | STOP |
---|
| 128 | ENDIF |
---|
| 129 | |
---|
| 130 | zciint = (zgdsup - za) / (zb - za) |
---|
| 131 | |
---|
| 132 | ! Computes the specific volume reference in situ temperature |
---|
| 133 | |
---|
| 134 | DO jk = 1, jpk |
---|
| 135 | zp = 0.e0 |
---|
[460] | 136 | zh = gdept_0(jk) |
---|
[3] | 137 | zt = 0.e0 |
---|
| 138 | zs = 35. |
---|
| 139 | zxk= zh * fsatg( zs, zt, zp ) |
---|
| 140 | zt = zt + 0.5 * zxk |
---|
| 141 | zq = zxk |
---|
| 142 | zp = zp + 0.5 * zh |
---|
| 143 | zxk= zh*fsatg( zs, zt, zp ) |
---|
| 144 | zt = zt + 0.29289322 * ( zxk - zq ) |
---|
| 145 | zq = 0.58578644 * zxk + 0.121320344 * zq |
---|
| 146 | zxk= zh * fsatg( zs, zt, zp ) |
---|
| 147 | zt = zt + 1.707106781 * ( zxk - zq ) |
---|
| 148 | zq = 3.414213562 * zxk - 4.121320344 * zq |
---|
| 149 | zp = zp + 0.5 * zh |
---|
| 150 | zxk= zh * fsatg( zs, zt, zp ) |
---|
| 151 | zwkx(jk) = zt + ( zxk - 2.0 * zq ) / 6.0 |
---|
| 152 | END DO |
---|
| 153 | |
---|
| 154 | ! In situ density (add the compression terms) |
---|
| 155 | |
---|
| 156 | DO jk = 1, jpk |
---|
| 157 | zt = zwkx(jk) |
---|
| 158 | zs = 35. |
---|
| 159 | ! square root salinity |
---|
| 160 | zsr = sqrt( abs( zs ) ) |
---|
| 161 | zwky(jk) = zsr |
---|
| 162 | ! compute density pure water at atm pressure |
---|
| 163 | zr1= ((((6.536332e-9*zt-1.120083e-6)*zt+1.001685e-4)*zt & |
---|
| 164 | -9.095290e-3)*zt+6.793952e-2)*zt+999.842594 |
---|
| 165 | ! seawater density atm pressure |
---|
| 166 | zr2= (((5.3875e-9*zt-8.2467e-7)*zt+7.6438e-5)*zt & |
---|
| 167 | -4.0899e-3)*zt+8.24493e-1 |
---|
| 168 | zr3= (-1.6546e-6*zt+1.0227e-4)*zt-5.72466e-3 |
---|
| 169 | zr4= 4.8314e-4 |
---|
| 170 | zwkz(jk)= (zr4*zs + zr3*zsr + zr2)*zs + zr1 |
---|
| 171 | END DO |
---|
| 172 | |
---|
| 173 | DO jk = 1, jpk |
---|
| 174 | zt = zwkx(jk) |
---|
| 175 | zs = 35. |
---|
| 176 | zsr= zwky(jk) |
---|
[460] | 177 | zh = gdept_0(jk) |
---|
[3] | 178 | |
---|
| 179 | ze = ( 9.1697e-11*zt+2.0816e-9 ) *zt-9.9348e-8 |
---|
| 180 | zbw= ( 5.2787e-9*zt-6.12293e-7 ) * zt+8.50935e-6 |
---|
| 181 | zb = zbw + ze * zs |
---|
| 182 | |
---|
| 183 | zd = 1.91075e-4 |
---|
| 184 | zc = (-1.6078e-6*zt-1.0981e-5)*zt+2.2838e-3 |
---|
| 185 | zaw= ((-5.77905e-7*zt+1.16092e-4)*zt+1.43713e-3)*zt+3.239908 |
---|
| 186 | za = ( zd*zsr + zc)*zs + zaw |
---|
| 187 | |
---|
| 188 | zb1= (-5.3009e-3*zt+1.6483e-1)*zt+7.944e-1 |
---|
| 189 | za1= ((-6.1670e-4*zt+1.09987e-1)*zt-6.03459)*zt+546.746 |
---|
| 190 | zkw= (((-5.155288e-4*zt+1.360477e-1)*zt-23.27105)*zt & |
---|
| 191 | +1484.206)*zt+196522.1 |
---|
| 192 | zk0= (zb1*zsr + za1)*zs + zkw |
---|
| 193 | ! evaluate pressure polynomial |
---|
| 194 | zwkz(jk) = zwkz(jk) / ( 1.0 - zh / ( zk0+zh*(za+zb*zh) ) ) |
---|
| 195 | END DO |
---|
| 196 | |
---|
| 197 | DO jk = 1, jpk |
---|
| 198 | rhosp(jk) = zwkz(jk) |
---|
| 199 | END DO |
---|
| 200 | ENDIF |
---|
| 201 | |
---|
| 202 | ! Computes the specific volume anomaly |
---|
| 203 | |
---|
| 204 | DO jk = 1, jpkm1 |
---|
| 205 | DO jj = 1, jpj |
---|
| 206 | DO ji = 1, jpi |
---|
| 207 | IF( tmask(ji,jj,jk) /= 0. ) THEN |
---|
| 208 | zsva(ji,jj,jk) = ( rau0*rhd(ji,jj,jk)+rau0 -rhosp(jk) ) / rhosp(jk) |
---|
| 209 | ELSE |
---|
| 210 | zsva(ji,jj,jk)=0. |
---|
| 211 | ENDIF |
---|
| 212 | END DO |
---|
| 213 | END DO |
---|
| 214 | END DO |
---|
| 215 | |
---|
| 216 | ! zfacto coefficient to cmg |
---|
| 217 | |
---|
| 218 | ! zfacto= 1. e+2 |
---|
| 219 | ! mg->cmg |
---|
| 220 | zfacto = 1.0 * 1.e2 |
---|
| 221 | |
---|
| 222 | ! Fisrt compute at depth ik=ihdsup |
---|
| 223 | |
---|
| 224 | ik = ihdsup |
---|
| 225 | DO jj = 1, jpj |
---|
| 226 | DO ji = 1, jpi |
---|
[460] | 227 | zhd = zfacto * zciint * e3t_0(ik) * zsva(ji,jj,ik) |
---|
[3] | 228 | hdy(ji,jj,ik) = zhd * tmask(ji,jj,ik) * tmask(ji,jj,ik-1) |
---|
| 229 | END DO |
---|
| 230 | END DO |
---|
| 231 | |
---|
| 232 | ! Then compute other terms except level jk=1 |
---|
| 233 | |
---|
| 234 | DO jk = ihdsup-1, 2, -1 |
---|
| 235 | DO jj = 1, jpj |
---|
| 236 | DO ji = 1, jpi |
---|
[460] | 237 | zhd = hdy(ji,jj,jk+1) + zfacto * e3t_0(jk) * zsva(ji,jj,jk) |
---|
[3] | 238 | hdy(ji,jj,jk) = zhd * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) |
---|
| 239 | END DO |
---|
| 240 | END DO |
---|
| 241 | END DO |
---|
| 242 | |
---|
| 243 | ! Then compute other the last layer term jk=1 |
---|
| 244 | |
---|
| 245 | ik = 1 |
---|
| 246 | DO jj = 1, jpj |
---|
| 247 | DO ji = 1, jpi |
---|
[460] | 248 | zhd = hdy(ji,jj,ik+1) + zfacto * e3t_0(ik) * zsva(ji,jj,ik) |
---|
[3] | 249 | hdy(ji,jj,ik) = zhd * tmask(ji,jj,ik) |
---|
| 250 | END DO |
---|
| 251 | END DO |
---|
| 252 | |
---|
| 253 | END SUBROUTINE dia_hdy |
---|
| 254 | |
---|
| 255 | #else |
---|
| 256 | !!---------------------------------------------------------------------- |
---|
| 257 | !! Default option : NO dynamic heigh diagnostics |
---|
| 258 | !!---------------------------------------------------------------------- |
---|
[32] | 259 | LOGICAL, PUBLIC, PARAMETER :: lk_diahdy = .FALSE. !: dynamical heigh flag |
---|
[3] | 260 | CONTAINS |
---|
| 261 | SUBROUTINE dia_hdy( kt ) ! Empty routine |
---|
[32] | 262 | WRITE(*,*) 'diahdy: You should not have seen this print! error?', kt |
---|
[3] | 263 | END SUBROUTINE dia_hdy |
---|
| 264 | #endif |
---|
| 265 | |
---|
| 266 | !!====================================================================== |
---|
| 267 | END MODULE diahdy |
---|