- Timestamp:
- 2011-09-27T14:33:01+02:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcvol.F90
r2797 r2865 10 10 !!---------------------------------------------------------------------- 11 11 #if defined key_obc && defined key_dynspg_flt 12 !!$ !!---------------------------------------------------------------------- 13 !!$ !! 'key_obc' AND unstructured open boundary conditions 14 !!$ !! 'key_dynspg_flt' filtered free surface 15 !!$ !!---------------------------------------------------------------------- 16 !!$ USE oce ! ocean dynamics and tracers 17 !!$ USE dom_oce ! ocean space and time domain 18 !!$ USE phycst ! physical constants 19 !!$ USE obc_oce ! ocean open boundary conditions 20 !!$ USE lib_mpp ! for mppsum 21 !!$ USE in_out_manager ! I/O manager 22 !!$ USE sbc_oce ! ocean surface boundary conditions 23 !!$ 24 !!$ IMPLICIT NONE 25 !!$ PRIVATE 26 !!$ 27 !!$ PUBLIC obc_vol ! routine called by dynspg_flt.h90 28 !!$ 29 !!$ !! * Substitutions 30 !!$# include "domzgr_substitute.h90" 31 !!$ !!---------------------------------------------------------------------- 32 !!$ !! NEMO/OPA 3.3 , NEMO Consortium (2010) 33 !!$ !! $Id$ 34 !!$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 35 !!$ !!---------------------------------------------------------------------- 36 !!$CONTAINS 37 !!$ 38 !!$ SUBROUTINE obc_vol( kt ) 39 !!$ !!---------------------------------------------------------------------- 40 !!$ !! *** ROUTINE obcvol *** 41 !!$ !! 42 !!$ !! ** Purpose : This routine is called in dynspg_flt to control 43 !!$ !! the volume of the system. A correction velocity is calculated 44 !!$ !! to correct the total transport through the unstructured OBC. 45 !!$ !! The total depth used is constant (H0) to be consistent with the 46 !!$ !! linear free surface coded in OPA 8.2 47 !!$ !! 48 !!$ !! ** Method : The correction velocity (zubtpecor here) is defined calculating 49 !!$ !! the total transport through all open boundaries (trans_obc) minus 50 !!$ !! the cumulate E-P flux (z_cflxemp) divided by the total lateral 51 !!$ !! surface (obcsurftot) of the unstructured boundary. 52 !!$ !! zubtpecor = [trans_obc - z_cflxemp ]*(1./obcsurftot) 53 !!$ !! with z_cflxemp => sum of (Evaporation minus Precipitation) 54 !!$ !! over all the domain in m3/s at each time step. 55 !!$ !! z_cflxemp < 0 when precipitation dominate 56 !!$ !! z_cflxemp > 0 when evaporation dominate 57 !!$ !! 58 !!$ !! There are 2 options (user's desiderata): 59 !!$ !! 1/ The volume changes according to E-P, this is the default 60 !!$ !! option. In this case the cumulate E-P flux are setting to 61 !!$ !! zero (z_cflxemp=0) to calculate the correction velocity. So 62 !!$ !! it will only balance the flux through open boundaries. 63 !!$ !! (set nn_volctl to 0 in tne namelist for this option) 64 !!$ !! 2/ The volume is constant even with E-P flux. In this case 65 !!$ !! the correction velocity must balance both the flux 66 !!$ !! through open boundaries and the ones through the free 67 !!$ !! surface. 68 !!$ !! (set nn_volctl to 1 in tne namelist for this option) 69 !!$ !!---------------------------------------------------------------------- 70 !!$ INTEGER, INTENT( in ) :: kt ! ocean time-step index 71 !!$ !! 72 !!$ INTEGER :: ji, jj, jk, jb, jgrd 73 !!$ INTEGER :: ii, ij 74 !!$ REAL(wp) :: zubtpecor, z_cflxemp, ztranst 75 !!$ !!----------------------------------------------------------------------------- 76 !!$ 77 !!$ IF( ln_vol ) THEN 78 !!$ 79 !!$ IF( kt == nit000 ) THEN 80 !!$ IF(lwp) WRITE(numout,*) 81 !!$ IF(lwp) WRITE(numout,*)'obc_vol : Correction of velocities along unstructured OBC' 82 !!$ IF(lwp) WRITE(numout,*)'~~~~~~~' 83 !!$ END IF 84 !!$ 85 !!$ ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain 86 !!$ ! ----------------------------------------------------------------------- 87 !!$ z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:) ) * obctmask(:,:) * e1t(:,:) * e2t(:,:) ) / rau0 88 !!$ IF( lk_mpp ) CALL mpp_sum( z_cflxemp ) ! sum over the global domain 89 !!$ 90 !!$ ! Transport through the unstructured open boundary 91 !!$ ! ------------------------------------------------ 92 !!$ zubtpecor = 0.e0 93 !!$ jgrd = 2 ! cumulate u component contribution first 94 !!$ DO jb = 1, nblenrim(jgrd) 95 !!$ DO jk = 1, jpkm1 96 !!$ ii = nbi(jb,jgrd) 97 !!$ ij = nbj(jb,jgrd) 98 !!$ zubtpecor = zubtpecor + flagu(jb) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 99 !!$ END DO 100 !!$ END DO 101 !!$ jgrd = 3 ! then add v component contribution 102 !!$ DO jb = 1, nblenrim(jgrd) 103 !!$ DO jk = 1, jpkm1 104 !!$ ii = nbi(jb,jgrd) 105 !!$ ij = nbj(jb,jgrd) 106 !!$ zubtpecor = zubtpecor + flagv(jb) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 107 !!$ END DO 108 !!$ END DO 109 !!$ IF( lk_mpp ) CALL mpp_sum( zubtpecor ) ! sum over the global domain 110 !!$ 111 !!$ ! The normal velocity correction 112 !!$ ! ------------------------------ 113 !!$ IF( nn_volctl==1 ) THEN ; zubtpecor = ( zubtpecor - z_cflxemp) / obcsurftot 114 !!$ ELSE ; zubtpecor = zubtpecor / obcsurftot 115 !!$ END IF 116 !!$ 117 !!$ ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation 118 !!$ ! ------------------------------------------------------------- 119 !!$ ztranst = 0.e0 120 !!$ jgrd = 2 ! correct u component 121 !!$ DO jb = 1, nblenrim(jgrd) 122 !!$ DO jk = 1, jpkm1 123 !!$ ii = nbi(jb,jgrd) 124 !!$ ij = nbj(jb,jgrd) 125 !!$ ua(ii,ij,jk) = ua(ii,ij,jk) - flagu(jb) * zubtpecor * umask(ii,ij,jk) 126 !!$ ztranst = ztranst + flagu(jb) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 127 !!$ END DO 128 !!$ END DO 129 !!$ jgrd = 3 ! correct v component 130 !!$ DO jb = 1, nblenrim(jgrd) 131 !!$ DO jk = 1, jpkm1 132 !!$ ii = nbi(jb,jgrd) 133 !!$ ij = nbj(jb,jgrd) 134 !!$ va(ii,ij,jk) = va(ii,ij,jk) -flagv(jb) * zubtpecor * vmask(ii,ij,jk) 135 !!$ ztranst = ztranst + flagv(jb) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 136 !!$ END DO 137 !!$ END DO 138 !!$ IF( lk_mpp ) CALL mpp_sum( ztranst ) ! sum over the global domain 139 !!$ 140 !!$ ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected 141 !!$ ! ------------------------------------------------------ 142 !!$ IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 143 !!$ IF(lwp) WRITE(numout,*) 144 !!$ IF(lwp) WRITE(numout,*)'obc_vol : time step :', kt 145 !!$ IF(lwp) WRITE(numout,*)'~~~~~~~ ' 146 !!$ IF(lwp) WRITE(numout,*)' cumulate flux EMP =', z_cflxemp , ' (m3/s)' 147 !!$ IF(lwp) WRITE(numout,*)' total lateral surface of OBC =', obcsurftot, '(m2)' 148 !!$ IF(lwp) WRITE(numout,*)' correction velocity zubtpecor =', zubtpecor , '(m/s)' 149 !!$ IF(lwp) WRITE(numout,*)' cumulated transport ztranst =', ztranst , '(m3/s)' 150 !!$ END IF 151 !!$ ! 152 !!$ END IF ! ln_vol 153 !!$ 154 !!$ END SUBROUTINE obc_vol 155 !!$ 156 !!$#else 12 !!---------------------------------------------------------------------- 13 !! 'key_obc' AND unstructured open boundary conditions 14 !! 'key_dynspg_flt' filtered free surface 15 !!---------------------------------------------------------------------- 16 USE oce ! ocean dynamics and tracers 17 USE dom_oce ! ocean space and time domain 18 USE phycst ! physical constants 19 USE obc_oce ! ocean open boundary conditions 20 USE lib_mpp ! for mppsum 21 USE in_out_manager ! I/O manager 22 USE sbc_oce ! ocean surface boundary conditions 23 24 IMPLICIT NONE 25 PRIVATE 26 27 PUBLIC obc_vol ! routine called by dynspg_flt.h90 28 29 !! * Substitutions 30 # include "domzgr_substitute.h90" 31 !!---------------------------------------------------------------------- 32 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 33 !! $Id$ 34 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 35 !!---------------------------------------------------------------------- 36 CONTAINS 37 38 SUBROUTINE obc_vol( kt ) 39 !!---------------------------------------------------------------------- 40 !! *** ROUTINE obcvol *** 41 !! 42 !! ** Purpose : This routine is called in dynspg_flt to control 43 !! the volume of the system. A correction velocity is calculated 44 !! to correct the total transport through the unstructured OBC. 45 !! The total depth used is constant (H0) to be consistent with the 46 !! linear free surface coded in OPA 8.2 47 !! 48 !! ** Method : The correction velocity (zubtpecor here) is defined calculating 49 !! the total transport through all open boundaries (trans_obc) minus 50 !! the cumulate E-P flux (z_cflxemp) divided by the total lateral 51 !! surface (obcsurftot) of the unstructured boundary. 52 !! zubtpecor = [trans_obc - z_cflxemp ]*(1./obcsurftot) 53 !! with z_cflxemp => sum of (Evaporation minus Precipitation) 54 !! over all the domain in m3/s at each time step. 55 !! z_cflxemp < 0 when precipitation dominate 56 !! z_cflxemp > 0 when evaporation dominate 57 !! 58 !! There are 2 options (user's desiderata): 59 !! 1/ The volume changes according to E-P, this is the default 60 !! option. In this case the cumulate E-P flux are setting to 61 !! zero (z_cflxemp=0) to calculate the correction velocity. So 62 !! it will only balance the flux through open boundaries. 63 !! (set nn_volctl to 0 in tne namelist for this option) 64 !! 2/ The volume is constant even with E-P flux. In this case 65 !! the correction velocity must balance both the flux 66 !! through open boundaries and the ones through the free 67 !! surface. 68 !! (set nn_volctl to 1 in tne namelist for this option) 69 !!---------------------------------------------------------------------- 70 INTEGER, INTENT( in ) :: kt ! ocean time-step index 71 !! 72 INTEGER :: ji, jj, jk, jb, jgrd 73 INTEGER :: ib_obc, ii, ij 74 REAL(wp) :: zubtpecor, z_cflxemp, ztranst 75 TYPE(OBC_INDEX), POINTER :: idx 76 !!----------------------------------------------------------------------------- 77 78 IF( ln_vol ) THEN 79 80 IF( kt == nit000 ) THEN 81 IF(lwp) WRITE(numout,*) 82 IF(lwp) WRITE(numout,*)'obc_vol : Correction of velocities along unstructured OBC' 83 IF(lwp) WRITE(numout,*)'~~~~~~~' 84 END IF 85 86 ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain 87 ! ----------------------------------------------------------------------- 88 z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:) ) * obctmask(:,:) * e1t(:,:) * e2t(:,:) ) / rau0 89 IF( lk_mpp ) CALL mpp_sum( z_cflxemp ) ! sum over the global domain 90 91 ! Transport through the unstructured open boundary 92 ! ------------------------------------------------ 93 zubtpecor = 0.e0 94 DO ib_obc = 1, nb_obc 95 idx => idx_obc(ib_obc) 96 97 jgrd = 2 ! cumulate u component contribution first 98 DO jb = 1, idx%nblenrim(jgrd) 99 DO jk = 1, jpkm1 100 ii = idx%nbi(jb,jgrd) 101 ij = idx%nbj(jb,jgrd) 102 zubtpecor = zubtpecor + idx%flagu(jb) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 103 END DO 104 END DO 105 jgrd = 3 ! then add v component contribution 106 DO jb = 1, idx%nblenrim(jgrd) 107 DO jk = 1, jpkm1 108 ii = idx%nbi(jb,jgrd) 109 ij = idx%nbj(jb,jgrd) 110 zubtpecor = zubtpecor + idx%flagv(jb) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 111 END DO 112 END DO 113 114 END DO 115 IF( lk_mpp ) CALL mpp_sum( zubtpecor ) ! sum over the global domain 116 117 ! The normal velocity correction 118 ! ------------------------------ 119 IF( nn_volctl==1 ) THEN ; zubtpecor = ( zubtpecor - z_cflxemp) / obcsurftot 120 ELSE ; zubtpecor = zubtpecor / obcsurftot 121 END IF 122 123 ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation 124 ! ------------------------------------------------------------- 125 ztranst = 0.e0 126 DO ib_obc = 1, nb_obc 127 idx => idx_obc(ib_obc) 128 129 jgrd = 2 ! correct u component 130 DO jb = 1, idx%nblenrim(jgrd) 131 DO jk = 1, jpkm1 132 ii = idx%nbi(jb,jgrd) 133 ij = idx%nbj(jb,jgrd) 134 ua(ii,ij,jk) = ua(ii,ij,jk) - idx%flagu(jb) * zubtpecor * umask(ii,ij,jk) 135 ztranst = ztranst + idx%flagu(jb) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 136 END DO 137 END DO 138 jgrd = 3 ! correct v component 139 DO jb = 1, idx%nblenrim(jgrd) 140 DO jk = 1, jpkm1 141 ii = idx%nbi(jb,jgrd) 142 ij = idx%nbj(jb,jgrd) 143 va(ii,ij,jk) = va(ii,ij,jk) -idx%flagv(jb) * zubtpecor * vmask(ii,ij,jk) 144 ztranst = ztranst + idx%flagv(jb) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 145 END DO 146 END DO 147 148 END DO 149 IF( lk_mpp ) CALL mpp_sum( ztranst ) ! sum over the global domain 150 151 ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected 152 ! ------------------------------------------------------ 153 IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 154 IF(lwp) WRITE(numout,*) 155 IF(lwp) WRITE(numout,*)'obc_vol : time step :', kt 156 IF(lwp) WRITE(numout,*)'~~~~~~~ ' 157 IF(lwp) WRITE(numout,*)' cumulate flux EMP =', z_cflxemp , ' (m3/s)' 158 IF(lwp) WRITE(numout,*)' total lateral surface of OBC =', obcsurftot, '(m2)' 159 IF(lwp) WRITE(numout,*)' correction velocity zubtpecor =', zubtpecor , '(m/s)' 160 IF(lwp) WRITE(numout,*)' cumulated transport ztranst =', ztranst , '(m3/s)' 161 END IF 162 ! 163 END IF ! ln_vol 164 165 END SUBROUTINE obc_vol 166 167 #else 157 168 !!---------------------------------------------------------------------- 158 169 !! Dummy module NO Unstruct Open Boundary Conditions
Note: See TracChangeset
for help on using the changeset viewer.
