[3611] | 1 | MODULE paresp |
---|
| 2 | !!====================================================================== |
---|
| 3 | !! *** MODULE paresp *** |
---|
| 4 | !! NEMOVAR : Weights for an energy-type scalar product |
---|
| 5 | !!====================================================================== |
---|
| 6 | |
---|
| 7 | !!---------------------------------------------------------------------- |
---|
| 8 | !! par_esp : Compute and store weights for an energy-type scalar product |
---|
| 9 | !!---------------------------------------------------------------------- |
---|
| 10 | !! * Modules used |
---|
| 11 | USE par_kind |
---|
| 12 | USE par_oce |
---|
| 13 | USE eosbn2 |
---|
| 14 | USE phycst |
---|
| 15 | USE zdf_oce |
---|
| 16 | USE dom_oce |
---|
| 17 | USE lib_mpp ! MPP library |
---|
| 18 | USE in_out_manager ! I/O stuff |
---|
| 19 | USE wrk_nemo ! Memory Allocation |
---|
| 20 | |
---|
| 21 | IMPLICIT NONE |
---|
| 22 | |
---|
| 23 | !! * Routine accessibility |
---|
| 24 | PRIVATE |
---|
| 25 | |
---|
| 26 | PUBLIC & |
---|
| 27 | & par_esp !: Compute and store energy weights |
---|
| 28 | |
---|
| 29 | !! * Share module variables |
---|
| 30 | REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC :: & |
---|
| 31 | & wesp_t, & !: Normalized energy weights for temperature |
---|
| 32 | & wesp_s !: Normalized energy weights for salinity |
---|
| 33 | REAL(wp), PUBLIC :: & |
---|
| 34 | & wesp_u, & !: Normalized energy weights for velocity |
---|
| 35 | & wesp_ssh, & !: Normalized energy weights for SSH |
---|
| 36 | & wesp_tau, & !: Normalized energy weights for windstress |
---|
| 37 | & wesp_q, & !: Normalized energy weights for heat flux |
---|
| 38 | & wesp_emp !: Normalized energy weights for EmP |
---|
| 39 | |
---|
| 40 | CONTAINS |
---|
| 41 | INTEGER FUNCTION par_esp_alloc() |
---|
| 42 | !!---------------------------------------------------------------------- |
---|
| 43 | !! *** FUNCTION exa_mpl_alloc *** |
---|
| 44 | !!---------------------------------------------------------------------- |
---|
| 45 | ALLOCATE( wesp_t(jpk) , wesp_s(jpk) , STAT= par_esp_alloc ) |
---|
| 46 | ! |
---|
| 47 | IF( lk_mpp ) CALL mpp_sum ( par_esp_alloc ) |
---|
| 48 | IF( par_esp_alloc /= 0 ) CALL ctl_warn('par_esp_alloc: failed to allocate arrays') |
---|
| 49 | ! |
---|
| 50 | END FUNCTION par_esp_alloc |
---|
| 51 | |
---|
| 52 | SUBROUTINE par_esp |
---|
| 53 | !!----------------------------------------------------------------------- |
---|
| 54 | !! |
---|
| 55 | !! *** ROUTINE par_esp *** |
---|
| 56 | !! |
---|
| 57 | !! ** Purpose : Compute and store weights for an energy-type scalar |
---|
| 58 | !! product. |
---|
| 59 | !! |
---|
| 60 | !! ** Method : |
---|
| 61 | !! The weighting coefficients are computed from an analytical |
---|
| 62 | !! linear density profile which is similar to the one used in |
---|
| 63 | !! default option in OPA direct. |
---|
| 64 | !! |
---|
| 65 | !! The weighting coefficients for the forcing fields are estimated |
---|
| 66 | !! from the surface boundary conditions. |
---|
| 67 | !! |
---|
| 68 | !! The weighting coefficients for the velocity variables are |
---|
| 69 | !! normalized to one. |
---|
| 70 | !! |
---|
| 71 | !! wesp_t(k) = ( g / N(k) )^2 x alpha^2 |
---|
| 72 | !! wesp_s(k) = ( g / N(k) )^2 x beta^2 |
---|
| 73 | !! wesp_u = 1 |
---|
| 74 | !! wesp_ssh = g |
---|
| 75 | !! wesp_tau = ( ( e3w(1) / ( avu x rho0 ) )^2 |
---|
| 76 | !! wesp_q = ( ( g / N(1) x alpha |
---|
| 77 | !! x e3w(1) / ( avT x rho0 x cp ) )^2 |
---|
| 78 | !! wesp_emp = ( ( g / N(1) x beta |
---|
| 79 | !! x e3w(1) x S(1) / avT )^2 |
---|
| 80 | !! |
---|
| 81 | !! where |
---|
| 82 | !! |
---|
| 83 | !! N(k) is the Brunt-Vaisala frequency (depth-dependent) |
---|
| 84 | !! alpha is the thermal expansion coefficient |
---|
| 85 | !! beta is the salinity expansion coefficient |
---|
| 86 | !! g is gravity |
---|
| 87 | !! avu is vertical eddy viscosity coefficent for dynamics |
---|
| 88 | !! avT is vertical eddy diffusivity coefficent for tracers |
---|
| 89 | !! e3w(1) is the vertical scale factor at top w-point |
---|
| 90 | !! rho0 is a reference density |
---|
| 91 | !! S(1) is a reference salinity at the surface |
---|
| 92 | !! |
---|
| 93 | !! The complete scalar product < ; > for the state vector |
---|
| 94 | !! |
---|
| 95 | !! dx = ( du , dv , dT , dS , deta , dq , demp , dtaux , dtauy ) |
---|
| 96 | !! |
---|
| 97 | !! involves the volume/area elements: |
---|
| 98 | !! |
---|
| 99 | !! < dx ; dx > = sum_i,j ( deta^2 * e1T * e2T * wesp_ssh |
---|
| 100 | !! + dq^2 * e1T * e2T * e3w(1) * wesp_q |
---|
| 101 | !! + demp^2 * e1T * e2T * e3w(1) * wesp_emp |
---|
| 102 | !! + dtaux^2 * e1u * e2u * e3uw(1) * wesp_tau |
---|
| 103 | !! + dtauy^2 * e1v * e2v * e3vw(1) * wesp_tau |
---|
| 104 | !! + sum_k ( du^2 * e1u * e2u * e3u * wesp_u |
---|
| 105 | !! + dv^2 * e1v * e2v * e3v * wesp_u |
---|
| 106 | !! + dT^2 * e1T * e2T * e3T * wesp_T |
---|
| 107 | !! + dS^2 * e1T * e2T * e3T * wesp_S ) ) |
---|
| 108 | !! |
---|
| 109 | !! ** Action : |
---|
| 110 | !! |
---|
| 111 | !! History : |
---|
| 112 | !! ! 97-06 (A. Weaver, J. Vialard) OPAVAR version |
---|
| 113 | !! ! 07-11 (A. Weaver) NEMOVAR version based on OPAVAR paresp.F |
---|
| 114 | !! ! 2011-07 (J. While) Adapted to deal with 2D grids |
---|
| 115 | !!----------------------------------------------------------------------- |
---|
| 116 | !! * Modules used |
---|
| 117 | |
---|
| 118 | !! * Arguments |
---|
| 119 | |
---|
| 120 | !! * Local declarations |
---|
| 121 | REAL(wp) :: zgrau0, zk, zt0, zs0, zrau0 |
---|
| 122 | REAL(wp), POINTER, DIMENSION(:) :: ztan, zsan |
---|
| 123 | REAL(wp), POINTER, DIMENSION(:) :: zrhan, zbn2an |
---|
| 124 | INTEGER :: jk, ierr |
---|
| 125 | !!---------------------------------------------------------------------- |
---|
| 126 | ! |
---|
| 127 | ierr = par_esp_alloc() |
---|
| 128 | CALL wrk_alloc( jpk, ztan, zsan, zrhan, zbn2an ) |
---|
| 129 | ! |
---|
| 130 | !-------------------------------------------------------------------- |
---|
| 131 | ! Local constant initialization |
---|
| 132 | !-------------------------------------------------------------------- |
---|
| 133 | |
---|
| 134 | zt0 = 19.0_wp ! Reference temperature |
---|
| 135 | zs0 = 35.0_wp ! Reference salinity |
---|
| 136 | zrau0 = 1025.0_wp ! Reference density |
---|
| 137 | |
---|
| 138 | zgrau0 = -1.0_wp * grav / zrau0 |
---|
| 139 | |
---|
| 140 | !-------------------------------------------------------------------- |
---|
| 141 | ! Initialize rho with an analytical profile |
---|
| 142 | !-------------------------------------------------------------------- |
---|
| 143 | |
---|
| 144 | ! Define analytical temperature profile |
---|
| 145 | |
---|
| 146 | DO jk = 1, jpk |
---|
| 147 | ztan(jk) = 7.5_wp * ( 1.0_wp - TANH( ( gdept_0(jk) - 80.0_wp ) & |
---|
| 148 | & / 30.0_wp ) ) & |
---|
| 149 | & + 10.0_wp * ( 6000.0_wp - gdept_0(jk) ) / 6000.0_wp |
---|
| 150 | ztan(jk) = ABS( ztan(jk) ) |
---|
| 151 | END DO |
---|
| 152 | #if defined key_pomme_r025 |
---|
| 153 | ztan(44) = 0.20 |
---|
| 154 | ztan(45) = 0.15 |
---|
| 155 | ztan(46) = 0.10 |
---|
| 156 | #endif |
---|
| 157 | |
---|
| 158 | ! Define analytical salinity profile |
---|
| 159 | |
---|
| 160 | DO jk = 1, jpk |
---|
| 161 | zsan(jk) = 35.50_wp |
---|
| 162 | END DO |
---|
| 163 | |
---|
| 164 | ! Define analytical density profile |
---|
| 165 | |
---|
| 166 | DO jk = 1, jpk |
---|
| 167 | zrhan(jk) = zrau0 * ( 1.0_wp - rn_alpha * ( ztan(jk) - zt0 ) & |
---|
| 168 | & + rn_beta * ( zsan(jk) - zs0 ) ) |
---|
| 169 | END DO |
---|
| 170 | |
---|
| 171 | !-------------------------------------------------------------------- |
---|
| 172 | ! Calculate N^2 at T and S points |
---|
| 173 | !-------------------------------------------------------------------- |
---|
| 174 | |
---|
| 175 | IF ( jpk > 2 ) THEN |
---|
| 176 | DO jk = 2, jpkm1 |
---|
| 177 | zbn2an(jk) = ( zgrau0 / 2.0_wp ) & |
---|
| 178 | & * ( ( ( zrhan(jk-1) - zrhan(jk ) ) / e3w_0(jk ) ) & |
---|
| 179 | & + ( ( zrhan(jk ) - zrhan(jk+1) ) / e3w_0(jk+1) ) ) |
---|
| 180 | END DO |
---|
| 181 | ELSE |
---|
| 182 | zbn2an(2) = 1._wp |
---|
| 183 | ENDIF |
---|
| 184 | |
---|
| 185 | zbn2an(1) = zbn2an(2) |
---|
| 186 | zbn2an(jpk) = zbn2an(jpkm1) |
---|
| 187 | |
---|
| 188 | !-------------------------------------------------------------------- |
---|
| 189 | ! Calculate energy-type weights |
---|
| 190 | !-------------------------------------------------------------------- |
---|
| 191 | |
---|
| 192 | DO jk = 1, jpk |
---|
| 193 | zk = ( grav * grav ) / zbn2an(jk) |
---|
| 194 | wesp_t(jk) = zk * rn_alpha * rn_alpha |
---|
| 195 | wesp_s(jk) = zk * rn_beta * rn_beta |
---|
| 196 | END DO |
---|
| 197 | |
---|
| 198 | wesp_u = 1.0_wp |
---|
| 199 | wesp_ssh = grav |
---|
| 200 | wesp_tau = ( e3w_0(1) / ( rn_avm0 * zrau0 ) )**2 |
---|
| 201 | wesp_q = wesp_t(1) * ( e3w_0(1) / ( rn_avt0 * zrau0 * rcp ) )**2 |
---|
| 202 | wesp_emp = wesp_s(1) * ( e3w_0(1) * zsan(1) / rn_avt0 )**2 |
---|
| 203 | |
---|
| 204 | !-------------------------------------------------------------------- |
---|
| 205 | ! Print |
---|
| 206 | !-------------------------------------------------------------------- |
---|
| 207 | |
---|
| 208 | IF (lwp) THEN |
---|
| 209 | WRITE(numout,*) |
---|
| 210 | WRITE(numout,*) ' par_esp: Analytical T, S, rho, N^2 profiles' |
---|
| 211 | WRITE(numout,*) ' -------' |
---|
| 212 | WRITE(numout,*) |
---|
| 213 | WRITE(numout,995) |
---|
| 214 | WRITE(numout,996) |
---|
| 215 | DO jk = 1, jpk |
---|
| 216 | IF(lwp)WRITE(numout,1007) jk, ztan(jk), zsan(jk), & |
---|
| 217 | & zrhan(jk), zbn2an(jk) |
---|
| 218 | END DO |
---|
| 219 | WRITE(numout,*) |
---|
| 220 | WRITE(numout,*) ' par_esp: Energy scalar product weights', & |
---|
| 221 | & ' for T and S' |
---|
| 222 | WRITE(numout,*) ' -------' |
---|
| 223 | WRITE(numout,*) |
---|
| 224 | WRITE(numout,998) |
---|
| 225 | WRITE(numout,999) |
---|
| 226 | DO jk = 1, jpk |
---|
| 227 | WRITE(numout,1006) jk, wesp_t(jk), wesp_s(jk) |
---|
| 228 | END DO |
---|
| 229 | WRITE(numout,*) |
---|
| 230 | WRITE(numout,*) ' par_esp: Energy scalar product weights', & |
---|
| 231 | & ' for u, SSH, tau, Q and EmP' |
---|
| 232 | WRITE(numout,*) ' -------' |
---|
| 233 | WRITE(numout,*) |
---|
| 234 | WRITE(numout,*) ' wesp_u = ', wesp_u |
---|
| 235 | WRITE(numout,*) ' wesp_ssh = ', wesp_ssh |
---|
| 236 | WRITE(numout,*) ' wesp_tau = ', wesp_tau |
---|
| 237 | WRITE(numout,*) ' wesp_q = ', wesp_q |
---|
| 238 | WRITE(numout,*) ' wesp_emp = ', wesp_emp |
---|
| 239 | |
---|
| 240 | 995 FORMAT(' level T S ', & |
---|
| 241 | & ' rho N^2 ') |
---|
| 242 | 996 FORMAT(' ----- -------- --------', & |
---|
| 243 | & ' ---------- ----------') |
---|
| 244 | 998 FORMAT(' level wesp_T wesp_S ') |
---|
| 245 | 999 FORMAT(' ----- ---------- ------------') |
---|
| 246 | 1005 FORMAT(4X,I2,3X,E12.5,1X,E12.5) |
---|
| 247 | 1006 FORMAT(4X,I2,3X,E12.5,1X,E12.5) |
---|
| 248 | 1007 FORMAT(4X,I2,6X,F10.5,5X,F10.5,5X,F10.5,5X,E12.5) |
---|
| 249 | CALL flush( numout ) |
---|
| 250 | |
---|
| 251 | ENDIF |
---|
| 252 | |
---|
| 253 | DO jk = 1, jpk |
---|
| 254 | IF ( zbn2an(jk) < 0.0_wp ) THEN |
---|
| 255 | IF(lwp)WRITE(numout,990) zbn2an(jk), jk |
---|
| 256 | 990 FORMAT(' paresp : Error unstable density profile;', & |
---|
| 257 | & ' zbn2an = ',E12.5,' at level ',I3) |
---|
| 258 | ! CALL ctl_stop( 'paresp; unstable density profile' ) |
---|
| 259 | IF(lwp)WRITE(numout,*)'WARNING' |
---|
| 260 | IF(lwp)WRITE(numout,*)'paresp; unstable density profile' |
---|
| 261 | ENDIF |
---|
| 262 | ENDDO |
---|
| 263 | CALL wrk_dealloc( jpk, ztan, zsan, zrhan, zbn2an ) |
---|
| 264 | |
---|
| 265 | END SUBROUTINE par_esp |
---|
| 266 | |
---|
| 267 | END MODULE paresp |
---|