Changeset 1125
- Timestamp:
- 2008-06-23T11:05:02+02:00 (16 years ago)
- Location:
- trunk/NEMO/OPA_SRC
- Files:
-
- 13 edited
-
BDY/bdy_oce.F90 (modified) (2 diffs)
-
BDY/bdy_par.F90 (modified) (1 diff)
-
BDY/bdydta.F90 (modified) (14 diffs)
-
BDY/bdydyn.F90 (modified) (6 diffs)
-
BDY/bdyini.F90 (modified) (7 diffs)
-
BDY/bdytides.F90 (modified) (8 diffs)
-
BDY/bdytra.F90 (modified) (2 diffs)
-
BDY/bdyvol.F90 (modified) (4 diffs)
-
DOM/domvvl.F90 (modified) (1 diff)
-
DYN/divcur.F90 (modified) (2 diffs)
-
DYN/dynnxt.F90 (modified) (1 diff)
-
DYN/wzvmod.F90 (modified) (2 diffs)
-
opa.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/BDY/bdy_oce.F90
r1058 r1125 4 4 !! Unstructured Open Boundary Cond. : define related variables 5 5 !!====================================================================== 6 #if defined key_bdy || defined key_bdy_tides 6 !! History : 1.0 ! 2001-05 (J. Chanut, A. Sellar) Original code 7 !! 3.0 ! 2008-04 (NEMO team) add in the reference version 7 8 !!---------------------------------------------------------------------- 8 !! 'key_bdy' : Unstructured Open Boundary Condition 9 #if defined key_bdy 9 10 !!---------------------------------------------------------------------- 10 !! History : 11 !! 9.0 01/05 (J. Chanut, A. Sellar) Original code 11 !! 'key_bdy' Unstructured Open Boundary Condition 12 12 !!---------------------------------------------------------------------- 13 !! * Modules used14 13 USE par_oce ! ocean parameters 15 14 USE bdy_par ! Unstructured boundary parameters … … 21 20 !! Namelist variables 22 21 !!---------------------------------------------------------------------- 23 24 CHARACTER(len=80) :: & 25 filbdy_mask, & !: Name of unstruct. bdy mask file 26 filbdy_data_T, & !: Name of unstruct. bdy data file at T points 27 filbdy_data_U, & !: Name of unstruct. bdy data file at U points 28 filbdy_data_V, & !: Name of unstruct. bdy data file at V points 29 filbdy_data_bt_T, & !: Name of unstruct. bdy data file at T points for barotropic variables 30 filbdy_data_bt_U, & !: Name of unstruct. bdy data file at U points for barotropic variables 31 filbdy_data_bt_V !: Name of unstruct. bdy data file at V points for barotropic variables 32 33 LOGICAL :: & 34 ln_bdy_clim =.false.!: if true, we assume that bdy data files contain 35 ! 1 time dump (-->bdy forcing will be constant) 36 ! or 12 months (-->bdy forcing will be cyclic) 37 LOGICAL :: & 38 ln_bdy_fla =.false. !: if true, flather boundary conditions 39 40 LOGICAL :: & 41 ln_bdy_vol =.false. !: if true, volume correction 42 43 LOGICAL :: & 44 ln_bdy_mask = .false. !: if true, read bdymask from file 45 46 INTEGER :: & 47 nb_rimwidth = 7, & !: boundary rim width 48 nbdy_dta = 1 !: = 0 use the initial state as bdy dta 49 !: = 1 read bdy data in netcdf file 50 INTEGER :: & 51 volbdy = 1 !: = 0 the total volume will have the variability 52 ! of the surface Flux E-P else (volbdy = 1) 53 ! the volume will be constant 54 ! = 1 the volume will be constant during all the 55 ! integration. 22 CHARACTER(len=80) :: filbdy_mask !: Name of unstruct. bdy mask file 23 CHARACTER(len=80) :: filbdy_data_T !: Name of unstruct. bdy data file at T points 24 CHARACTER(len=80) :: filbdy_data_U !: Name of unstruct. bdy data file at U points 25 CHARACTER(len=80) :: filbdy_data_V !: Name of unstruct. bdy data file at V points 26 CHARACTER(len=80) :: filbdy_data_bt_T !: Name of unstruct. bdy data file at T points for barotropic variables 27 CHARACTER(len=80) :: filbdy_data_bt_U !: Name of unstruct. bdy data file at U points for barotropic variables 28 CHARACTER(len=80) :: filbdy_data_bt_V !: Name of unstruct. bdy data file at V points for barotropic variables 29 ! 30 LOGICAL :: ln_bdy_tides = .false. !: =T apply tidal harmonic forcing along open boundaries 31 LOGICAL :: ln_bdy_vol = .false. !: =T volume correction 32 LOGICAL :: ln_bdy_mask = .false. !: =T read bdymask from file 33 LOGICAL :: ln_bdy_clim = .false. !: if true, we assume that bdy data files contain 34 ! ! 1 time dump (-->bdy forcing will be constant) 35 ! ! or 12 months (-->bdy forcing will be cyclic) 36 LOGICAL :: ln_bdy_dyn_fla = .false. !: =T Flather boundary conditions on barotropic velocities 37 LOGICAL :: ln_bdy_dyn_frs = .false. !: =T FRS boundary conditions on velocities 38 LOGICAL :: ln_bdy_tra_frs = .false. !: =T FRS boundary conditions on tracers (T and S) 39 LOGICAL :: ln_bdy_ice_frs = .false. !: =T FRS boundary conditions on seaice (leads fraction, ice depth, snow depth) 40 ! 41 INTEGER :: nb_rimwidth = 7 !: boundary rim width 42 INTEGER :: nbdy_dta = 1 !: = 0 use the initial state as bdy dta or = 1 read it in a NetCDF file 43 INTEGER :: volbdy = 1 !: = 0 the total volume will have the variability of the surface Flux E-P 44 ! ! = 1 the volume will be constant during all the integration. 56 45 57 46 !!---------------------------------------------------------------------- 58 47 !! Global variables 59 48 !!---------------------------------------------------------------------- 60 61 REAL(wp), DIMENSION(jpi, jpj) :: & !: 62 bdytmask, & !: Mask defining computational domain at T-points 63 bdyumask, & !: Mask defining computational domain at U-points 64 bdyvmask !: Mask defining computational domain at V-points 49 REAL(wp), DIMENSION(jpi,jpj) :: bdytmask !: Mask defining computational domain at T-points 50 REAL(wp), DIMENSION(jpi,jpj) :: bdyumask !: Mask defining computational domain at U-points 51 REAL(wp), DIMENSION(jpi,jpj) :: bdyvmask !: Mask defining computational domain at V-points 65 52 66 53 !!---------------------------------------------------------------------- 67 54 !! Unstructured open boundary data variables 68 55 !!---------------------------------------------------------------------- 56 INTEGER, DIMENSION(jpbgrd) :: nblen !: Size of bdy data on a proc for each grid type 57 INTEGER, DIMENSION(jpbgrd) :: nblenrim !: Size of bdy data on a proc for first rim ind 58 INTEGER, DIMENSION(jpbgrd) :: nblendta !: Size of bdy data in file 69 59 70 INTEGER, DIMENSION(jpbgrd) :: & 71 nblen = 0, & !: Size of bdy data on a proc for each grid type 72 nblenrim = 0, & !: Size of bdy data on a proc for first rim ind 73 nblendta = 0 !: Size of bdy data in file 60 INTEGER, DIMENSION(jpbdim,jpbgrd) :: nbi, nbj !: i and j indices of bdy dta 61 INTEGER, DIMENSION(jpbdim,jpbgrd) :: nbr !: Discrete distance from rim points 62 INTEGER, DIMENSION(jpbdim,jpbgrd) :: nbmap !: Indices of data in file for data in memory 63 64 REAL(wp) :: bdysurftot !: Lateral surface of unstructured open boundary 74 65 75 INTEGER, DIMENSION(jpbdim, jpbgrd) :: & 76 nbi, nbj, & !: i and j indices of bdy dta 77 nbr, & !: Discrete distance from rim points 78 nbmap !: Indices of data in file for data in memory 79 80 REAL(wp) :: & 81 bdysurftot !: Lateral surface of unstructured open boundary 66 REAL(wp), DIMENSION(jpbdim) :: flagu, flagv !: Flag for normal velocity compnt for velocity components 67 REAL(wp), DIMENSION(jpbdim,jpbgrd) :: nbw !: Rim weights of bdy data 82 68 83 REAL(wp), DIMENSION(jpbdim) :: & 84 flagu, & !: Flag for normal velocity compnt for u-velocity 85 flagv !: Flag for normal velocity compnt for v-velocity 69 REAL(wp), DIMENSION(jpbdim) :: sshbdy !: Now clim of bdy sea surface height (Flather) 70 REAL(wp), DIMENSION(jpbdim) :: ubtbdy, vbtbdy !: Now clim of bdy barotropic velocity components 71 REAL(wp), DIMENSION(jpbdim,jpk) :: tbdy , sbdy !: Now clim of bdy temperature and salinity 72 REAL(wp), DIMENSION(jpbdim,jpk) :: ubdy , vbdy !: Now clim of bdy velocity components 73 REAL(wp), DIMENSION(jpbdim) :: sshtide !: Tidal boundary array : SSH 74 REAL(wp), DIMENSION(jpbdim) :: utide, vtide !: Tidal boundary array : U and V 86 75 87 REAL(wp), DIMENSION(jpbdim, jpbgrd) :: &88 nbw !: Rim weights of bdy data89 90 REAL(wp), DIMENSION(jpbdim) :: &91 sshbdy, & !: Now clim of bdy sea surface height (Flather)92 ubtbdy, vbtbdy !: Now clim of bdy barotropic velocity components93 94 REAL(wp), DIMENSION(jpbdim,jpk) :: &95 tbdy, sbdy, & !: Now clim of bdy temperature and salinity96 ubdy, vbdy !: Now clim of bdy velocity components97 98 #if defined key_bdy_tides99 REAL(wp), DIMENSION(jpbdim) :: &100 sshtide, & !: Tidal boundary array : SSH101 utide, & !: Tidal boundary array : U102 vtide !: Tidal boundary array : V103 #endif104 105 76 #else 106 77 !!---------------------------------------------------------------------- 107 !! D efault option : Empty module78 !! Dummy module NO Unstructured Open Boundary Condition 108 79 !!---------------------------------------------------------------------- 109 80 #endif 110 81 82 !!---------------------------------------------------------------------- 83 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 84 !! $Id: $ 85 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 111 86 !!====================================================================== 112 87 END MODULE bdy_oce -
trunk/NEMO/OPA_SRC/BDY/bdy_par.F90
r1058 r1125 1 1 MODULE bdy_par 2 2 !!====================================================================== 3 !! *** MODULE bdy_par ***3 !! *** MODULE bdy_par *** 4 4 !! Unstructured Open Boundary Cond. : define related parameters 5 5 !!====================================================================== 6 #if defined key_bdy || defined key_bdy_tides 6 !! History : 1.0 ! 2005-01 (J. Chanut, A. Sellar) Original code 7 !! 3.0 ! 2008-04 (NEMO team) add in the reference version 8 !!---------------------------------------------------------------------- 9 #if defined key_bdy 7 10 !!---------------------------------------------------------------------- 8 11 !! 'key_bdy' : Unstructured Open Boundary Condition 9 12 !!---------------------------------------------------------------------- 10 !! history :11 !! 9.0 01/05 (J. Chanut, A. Sellar) Original code12 !!----------------------------------------------------------------------13 !! * Modules used14 USE par_oce ! ocean parameters15 13 16 14 IMPLICIT NONE 17 15 PUBLIC 18 16 19 #if defined key_bdy 20 LOGICAL, PUBLIC, PARAMETER :: lk_bdy = .TRUE. !: Unstructured Ocean 21 !Boundary Condition flag 17 18 LOGICAL, PUBLIC, PARAMETER :: lk_bdy = .TRUE. !: Unstructured Ocean Boundary Condition flag 19 INTEGER, PUBLIC, PARAMETER :: jpbdta = 5000 !: Max length of bdy field in file 20 INTEGER, PUBLIC, PARAMETER :: jpbdim = 5000 !: Max length of bdy field on a processor 21 INTEGER, PUBLIC, PARAMETER :: jpbtime = 1000 !: Max number of time dumps per file 22 INTEGER, PUBLIC, PARAMETER :: jpbgrd = 3 !: Number of horizontal grid types used (T, u, v, f) 22 23 #else 23 ! tides only at boundaries 24 LOGICAL, PUBLIC, PARAMETER :: lk_bdy = .FALSE. !: Unstructured Ocean 25 !Boundary Condition flag 24 !!---------------------------------------------------------------------- 25 !! Default option : NO Unstructured open boundary condition 26 !!---------------------------------------------------------------------- 27 LOGICAL, PUBLIC, PARAMETER :: lk_bdy = .FALSE. !: Unstructured Ocean Boundary Condition flag 26 28 #endif 27 29 28 30 !!---------------------------------------------------------------------- 29 !! Unstructured open boundary parameters 30 !!---------------------------------------------------------------------- 31 32 INTEGER, PARAMETER :: & 33 jpbdta = 5000, & !: Max length of bdy field in file 34 jpbdim = 5000, & !: Max length of bdy field on a processor 35 jpbtime = 1000, & !: Max number of time dumps per file 36 jpbgrd = 3 !: Number of horizontal grid types used 37 ! (T, u, v, f) 38 #else 39 !!---------------------------------------------------------------------- 40 !! Default option : NO open boundary condition 41 !!---------------------------------------------------------------------- 42 LOGICAL, PUBLIC, PARAMETER :: lk_bdy = .FALSE.!: Unstructured Ocean 43 !Boundary Condition flag 44 #endif 45 31 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 32 !! $Id: $ 33 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 46 34 !!====================================================================== 47 35 END MODULE bdy_par -
trunk/NEMO/OPA_SRC/BDY/bdydta.F90
r1058 r1125 1 1 MODULE bdydta 2 !!====================================================================== ========3 !! 2 !!====================================================================== 3 !! *** MODULE bdydta *** 4 4 !! Open boundary data : read the data for the unstructured open boundaries. 5 !!============================================================================== 6 #if defined key_bdy || defined key_bdy_tides 7 !!------------------------------------------------------------------------------ 8 !! 'key_bdy' : Unstructured Open Boundary Conditions 9 !!------------------------------------------------------------------------------ 10 !! bdy_dta : read u, v, t, s data along each open boundary 11 !!------------------------------------------------------------------------------ 12 !! * Modules used 5 !!====================================================================== 6 !! History : 1.0 ! 2005-01 (J. Chanut, A. Sellar) Original code 7 !! - ! 2007-01 (D. Storkey) Update to use IOM module 8 !! - ! 2007-07 (D. Storkey) add bdy_dta_bt 9 !! 3.0 ! 2008-04 (NEMO team) add in the reference version 10 !!---------------------------------------------------------------------- 11 #if defined key_bdy 12 !!---------------------------------------------------------------------- 13 !! 'key_bdy' Unstructured Open Boundary Conditions 14 !!---------------------------------------------------------------------- 15 !! bdy_dta : read u, v, t, s data along open boundaries 16 !! bdy_dta_bt : read depth-mean velocities and elevation along open 17 !! boundaries 18 !!---------------------------------------------------------------------- 13 19 USE oce ! ocean dynamics and tracers 14 20 USE dom_oce ! ocean space and time domain … … 24 30 PRIVATE 25 31 26 !! * Accessibility 27 PUBLIC bdy_dta ! routines called by step.F90 28 PUBLIC bdy_dta_bt 32 PUBLIC bdy_dta ! routines called by step.F90 33 PUBLIC bdy_dta_bt 34 35 INTEGER :: numbdyt, numbdyu, numbdyv !: logical units for T-, U-, & V-points data file, resp. 36 INTEGER :: ntimes_bdy !: exact number of time dumps in data files 37 INTEGER :: nbdy_b, nbdy_a !: record of bdy data file for before and after model time step 38 INTEGER :: numbdyt_bt, numbdyu_bt, numbdyv_bt !: logical unit for T-, U- & V-points data file, resp. 39 INTEGER :: ntimes_bdy_bt !: exact number of time dumps in data files 40 INTEGER :: nbdy_b_bt, nbdy_a_bt !: record of bdy data file for before and after model time step 41 42 INTEGER, DIMENSION (jpbtime) :: istep, istep_bt !: time array in seconds in each data file 43 44 REAL(wp) :: zoffset !: time offset between time origin in file & start time of model run 45 46 REAL(wp), DIMENSION(jpbdim,jpk,2) :: tbdydta, sbdydta !: time interpolated values of T and S bdy data 47 REAL(wp), DIMENSION(jpbdim,jpk,2) :: ubdydta, vbdydta !: time interpolated values of U and V bdy data 48 REAL(wp), DIMENSION(jpbdim,2) :: ubtbdydta, vbtbdydta !: Arrays used for time interpolation of bdy data 49 REAL(wp), DIMENSION(jpbdim,2) :: sshbdydta !: bdy data of ssh 50 51 !!---------------------------------------------------------------------- 52 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 53 !! $Id: $ 54 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 55 !!---------------------------------------------------------------------- 56 57 CONTAINS 58 59 SUBROUTINE bdy_dta( kt ) 60 !!---------------------------------------------------------------------- 61 !! *** SUBROUTINE bdy_dta *** 62 !! 63 !! ** Purpose : Read unstructured boundary data for FRS condition. 64 !! 65 !! ** Method : At the first timestep, read in boundary data for two 66 !! times from the file and time-interpolate. At other 67 !! timesteps, check to see if we need another time from 68 !! the file. If so read it in. Time interpolate. 69 !!---------------------------------------------------------------------- 70 INTEGER, INTENT( in ) :: kt ! ocean time-step index (for timesplitting option, otherwise zero) 71 !! 72 CHARACTER(LEN=80), DIMENSION(3) :: clfile ! names of input files 73 CHARACTER(LEN=70 ) :: clunits ! units attribute of time coordinate 74 LOGICAL :: lect ! flag for reading 75 INTEGER :: it, ib, ik, igrd ! dummy loop indices 76 INTEGER :: igrd_start, igrd_end ! start and end of loops on igrd 77 INTEGER :: idvar ! netcdf var ID 78 INTEGER :: iman, i15, imois ! Time variables for monthly clim forcing 79 INTEGER :: ntimes_bdyt, ntimes_bdyu, ntimes_bdyv 80 INTEGER :: itimer, totime 81 INTEGER :: ii, ij ! array addresses 82 INTEGER :: ipi, ipj, ipk, inum ! temporary integers (NetCDF read) 83 INTEGER :: iyear0, imonth0, iday0 84 INTEGER :: ihours0, iminutes0, isec0 85 INTEGER :: iyear, imonth, iday, isecs 86 INTEGER, DIMENSION(jpbtime) :: istept, istepu, istepv ! time arrays from data files 87 REAL(wp) :: dayfrac, zxy, zoffsett 88 REAL(wp) :: zoffsetu, zoffsetv 89 REAL(wp) :: dayjul0, zdayjulini 90 REAL(wp), DIMENSION(jpbtime) :: zstepr ! REAL time array from data files 91 REAL(wp), DIMENSION(jpbdta,1,jpk) :: zdta ! temporary array for data fields 92 !!--------------------------------------------------------------------------- 93 94 IF( ln_bdy_dyn_frs .OR. ln_bdy_tra_frs ) THEN ! If these are both false then this routine 95 ! does nothing. 96 97 ! -------------------- ! 98 ! Initialization ! 99 ! -------------------- ! 100 101 lect = .false. ! If true, read a time record 102 103 ! Some time variables for monthly climatological forcing: 104 ! ******************************************************* 105 !!gm here use directely daymod variables 106 107 iman = INT( raamo ) ! Number of months in a year 108 109 i15 = INT( 2*FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) ) 110 ! i15=0 if the current day is in the first half of the month, else i15=1 111 112 imois = nmonth + i15 - 1 ! imois is the first month record 113 IF( imois == 0 ) imois = iman 114 115 ! Time variable for non-climatological forcing: 116 ! ********************************************* 117 itimer = (kt-nit000+1)*rdt ! current time in seconds for interpolation 118 119 120 ! !-------------------! 121 IF( kt == nit000 ) THEN ! First call only ! 122 ! !-------------------! 123 istep(:) = 0 124 nbdy_b = 0 125 nbdy_a = 0 126 127 ! Get time information from bdy data file 128 ! *************************************** 129 130 IF(lwp) WRITE(numout,*) 131 IF(lwp) WRITE(numout,*) 'bdy_dta : Initialize unstructured boundary data' 132 IF(lwp) WRITE(numout,*) '~~~~~~~' 133 134 IF ( nbdy_dta == 0 ) THEN 135 ! 136 IF(lwp) WRITE(numout,*) ' Bdy data are taken from initial conditions' 137 ! 138 ELSEIF (nbdy_dta == 1) THEN 139 ! 140 IF(lwp) WRITE(numout,*) ' Bdy data are read in netcdf files' 141 ! 142 dayfrac = adatrj - FLOAT( itimer ) / 86400. ! day fraction at time step kt-1 143 dayfrac = dayfrac - INT ( dayfrac ) ! 144 totime = ( nitend - nit000 + 1 ) * rdt ! Total time of the run to verify that all the 145 ! ! necessary time dumps in file are included 146 ! 147 clfile(1) = filbdy_data_T 148 clfile(2) = filbdy_data_U 149 clfile(3) = filbdy_data_V 150 ! 151 ! how many files are we to read in? 152 igrd_start = 1 153 igrd_end = 3 154 IF(.NOT. ln_bdy_tra_frs .AND. .NOT. ln_bdy_ice_frs) THEN 155 ! No T-grid file. 156 igrd_start = 2 157 ELSEIF ( .NOT. ln_bdy_dyn_frs ) THEN 158 ! No U-grid or V-grid file. 159 igrd_end = 1 160 ENDIF 161 162 DO igrd = igrd_start, igrd_end ! loop over T, U & V grid ! 163 ! !---------------------------! 164 CALL iom_open( clfile(igrd), inum ) 165 CALL iom_gettime( inum, zstepr, kntime=ntimes_bdy, cdunits=clunits ) 166 167 SELECT CASE( igrd ) 168 CASE (1) 169 numbdyt = inum 170 CASE (2) 171 numbdyu = inum 172 CASE (3) 173 numbdyv = inum 174 END SELECT 175 176 ! Calculate time offset 177 READ(clunits,7000) iyear0, imonth0, iday0, ihours0, iminutes0, isec0 178 ! Convert time origin in file to julian days 179 isec0 = isec0 + ihours0*60.*60. + iminutes0*60. 180 CALL ymds2ju(iyear0, imonth0, iday0, real(isec0), dayjul0) 181 ! Compute model initialization time 182 iyear = ndastp / 10000 183 imonth = ( ndastp - iyear * 10000 ) / 100 184 iday = ndastp - iyear * 10000 - imonth * 100 185 isecs = dayfrac * 86400 186 CALL ymds2ju(iyear, imonth, iday, real(isecs) , zdayjulini) 187 ! offset from initialization date: 188 zoffset = (dayjul0-zdayjulini)*86400 189 ! 190 7000 FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 191 192 !! TO BE DONE... Check consistency between calendar from file 193 !! (available optionally from iom_gettime) and calendar in model 194 !! when calendar in model available outside of IOIPSL. 195 196 IF(lwp) WRITE(numout,*) 'number of times: ',ntimes_bdy 197 IF(lwp) WRITE(numout,*) 'offset: ',zoffset 198 IF(lwp) WRITE(numout,*) 'totime: ',totime 199 IF(lwp) WRITE(numout,*) 'zstepr: ',zstepr 200 201 ! Check that there are not too many times in the file. 202 IF( ntimes_bdy > jpbtime ) THEN 203 WRITE(ctmp1,*) 'Check file: ', clfile(igrd), 'jpbtime= ', jpbtime, ' ntimes_bdy= ', ntimes_bdy 204 CALL ctl_stop( 'Number of time dumps in files exceed jpbtime parameter', ctmp1 ) 205 ENDIF 206 207 ! Check that time array increases: 208 209 it = 1 210 DO WHILE( zstepr(it+1) > zstepr(it) .AND. it /= ntimes_bdy - 1 ) 211 it = it + 1 212 END DO 213 214 IF( it.NE.ntimes_bdy-1 .AND. ntimes_bdy > 1 ) THEN 215 WRITE(ctmp1,*) 'Check file: ', clfile(igrd) 216 CALL ctl_stop( 'Time array in unstructured boundary data files', & 217 & 'does not continuously increase.' , ctmp1 ) 218 ENDIF 219 ! 220 ! Check that times in file span model run time: 221 IF( zstepr(1) + zoffset > 0 ) THEN 222 WRITE(ctmp1,*) 'Check file: ', clfile(igrd) 223 CALL ctl_stop( 'First time dump in bdy file is after model initial time', ctmp1 ) 224 END IF 225 IF( zstepr(ntimes_bdy) + zoffset < totime ) THEN 226 WRITE(ctmp1,*) 'Check file: ', clfile(igrd) 227 CALL ctl_stop( 'Last time dump in bdy file is before model final time', ctmp1 ) 228 END IF 229 ! 230 IF ( igrd == 1 ) THEN 231 ntimes_bdyt = ntimes_bdy 232 zoffsett = zoffset 233 istept(:) = INT( zstepr(:) + zoffset ) 234 ELSEIF(igrd == 2 ) THEN 235 ntimes_bdyu = ntimes_bdy 236 zoffsetu = zoffset 237 istepu(:) = INT( zstepr(:) + zoffset ) 238 ELSEIF(igrd == 3 ) THEN 239 ntimes_bdyv = ntimes_bdy 240 zoffsetv = zoffset 241 istepv(:) = INT( zstepr(:) + zoffset ) 242 ENDIF 243 ! 244 END DO ! end loop over T, U & V grid 245 246 IF (igrd_start == 1 .and. igrd_end == 3) THEN 247 ! Only test differences if we are reading in 3 files 248 ! Verify time consistency between files 249 IF( ntimes_bdyu /= ntimes_bdyt .OR. ntimes_bdyv /= ntimes_bdyt ) THEN 250 CALL ctl_stop( 'Bdy data files must have the same number of time dumps', & 251 & 'Multiple time frequencies not implemented yet' ) 252 ENDIF 253 ntimes_bdy = ntimes_bdyt 254 ! 255 IF( zoffsetu /= zoffsett .OR. zoffsetv /= zoffsett ) THEN 256 CALL ctl_stop( 'Bdy data files must have the same time origin', & 257 & 'Multiple time frequencies not implemented yet' ) 258 ENDIF 259 zoffset = zoffsett 260 ENDIF 261 262 IF( igrd_start == 1 ) THEN 263 istep(:) = istept(:) 264 ELSE 265 istep(:) = istepu(:) 266 ENDIF 267 268 ! Check number of time dumps: 269 IF( ntimes_bdy == 1 .AND. .NOT. ln_bdy_clim ) THEN 270 CALL ctl_stop( 'There is only one time dump in data files', & 271 & 'Choose ln_bdy_clim=.true. in namelist for constant bdy forcing.' ) 272 ENDIF 273 274 IF( ln_bdy_clim ) THEN 275 IF( ntimes_bdy /= 1 .AND. ntimes_bdy /= 12 ) THEN 276 CALL ctl_stop( 'For climatological boundary forcing (ln_bdy_clim=.true.),', & 277 & 'bdy data files must contain 1 or 12 time dumps.' ) 278 ELSEIF( ntimes_bdy == 1 ) THEN 279 IF(lwp) WRITE(numout,*) 280 IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files' 281 ELSEIF( ntimes_bdy == 12 ) THEN 282 IF(lwp) WRITE(numout,*) 283 IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files' 284 ENDIF 285 ENDIF 286 287 ! Find index of first record to read (before first model time). 288 it = 1 289 DO WHILE( istep(it+1) <= 0 .AND. it <= ntimes_bdy - 1 ) 290 it = it + 1 291 END DO 292 nbdy_b = it 293 ! 294 WRITE(numout,*) 'Time offset is ',zoffset 295 WRITE(numout,*) 'First record to read is ',nbdy_b 296 297 ENDIF ! endif (nbdy_dta == 1) 298 299 300 ! 1.2 Read first record in file if necessary (ie if nbdy_dta == 1) 301 ! ***************************************************************** 302 303 IF( nbdy_dta == 0) THEN ! boundary data arrays are filled with initial conditions 304 ! 305 IF (ln_bdy_tra_frs) THEN 306 igrd = 1 ! T-points data 307 DO ib = 1, nblen(igrd) 308 ii = nbi(ib,igrd) 309 ij = nbj(ib,igrd) 310 DO ik = 1, jpkm1 311 tbdy(ib,ik) = tn(ii, ij, ik) 312 sbdy(ib,ik) = sn(ii, ij, ik) 313 ENDDO 314 END DO 315 ENDIF 316 317 IF(ln_bdy_dyn_frs) THEN 318 igrd = 2 ! U-points data 319 DO ib = 1, nblen(igrd) 320 ii = nbi(ib,igrd) 321 ij = nbj(ib,igrd) 322 DO ik = 1, jpkm1 323 ubdy(ib,ik) = un(ii, ij, ik) 324 ENDDO 325 END DO 326 327 igrd = 3 ! V-points data 328 DO ib = 1, nblen(igrd) 329 ii = nbi(ib,igrd) 330 ij = nbj(ib,igrd) 331 DO ik = 1, jpkm1 332 vbdy(ib,ik) = vn(ii, ij, ik) 333 ENDDO 334 END DO 335 ENDIF 336 ! 337 ELSEIF( nbdy_dta == 1 ) THEN ! Set first record in the climatological case: 338 ! 339 IF( ln_bdy_clim .AND. ntimes_bdy == 1 ) THEN 340 nbdy_a = 1 341 ELSEIF( ln_bdy_clim .AND. ntimes_bdy == iman ) THEN 342 nbdy_b = 0 343 nbdy_a = imois 344 ELSE 345 nbdy_a = nbdy_b 346 ENDIF 29 347 30 !! * local modules variables 31 INTEGER :: & 32 bdynumt, & ! logical unit for T-points data file 33 bdynumu, & ! logical unit for U-points data file 34 bdynumv, & ! logical unit for V-points data file 35 kbdy, & ! Exact number of time dumps in data files 36 nbdy1, & ! Time record in bdy data file BEFORE the model current time 37 nbdy2 ! Time record in bdy data file AFTER the model current time 38 39 !! * local modules variables for barotropic variables 40 INTEGER :: & 41 bdynumt_bt, & ! logical unit for T-points data file 42 bdynumu_bt, & ! logical unit for U-points data file 43 bdynumv_bt, & ! logical unit for V-points data file 44 kbdy_bt, & ! Exact number of time dumps in data files 45 nbdy1_bt, & ! Time record in bdy data file BEFORE the model current time 46 nbdy2_bt ! Time record in bdy data file AFTER the model current time 47 48 INTEGER, DIMENSION (jpbtime) :: & 49 istep, istep_bt ! time array in seconds in each data file 50 51 REAL(wp) :: & 52 offset ! time offset between time origin in file and start time of model run 53 54 REAL(wp), DIMENSION(jpbdim,jpk,2) :: & 55 tbdydta, sbdydta, & !: Arrays used for time interpolation of bdy data 56 ubdydta, vbdydta 57 58 REAL(wp), DIMENSION(jpbdim,2) :: & 59 ubtbdydta, vbtbdydta, & !: Arrays used for time interpolation of bdy data 60 sshbdydta !: bdy data of ssh 61 62 !!--------------------------------------------------------------------------------- 63 !! OPA 9.0 , LODYC-IPSL (2003) 64 !!--------------------------------------------------------------------------------- 65 66 CONTAINS 67 68 SUBROUTINE bdy_dta( kt ) 348 ! Read first record: 349 ipj = 1 350 ipk = jpk 351 igrd = 1 352 ipi = nblendta(igrd) 353 354 IF(ln_bdy_tra_frs) THEN 355 igrd = 1 ! Temperature 356 IF( nblendta(igrd) <= 0 ) THEN 357 idvar = iom_varid( numbdyt, 'votemper' ) 358 nblendta(igrd) = iom_file(numbdyt)%dimsz(1,idvar) 359 ENDIF 360 WRITE(numout,*) 'Dim size for votemper is ', nblendta(igrd) 361 ipi = nblendta(igrd) 362 CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 363 364 DO ib = 1, nblen(igrd) 365 DO ik = 1, jpkm1 366 tbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 367 END DO 368 END DO 369 ! 370 igrd = 1 ! salinity 371 IF( nblendta(igrd) .le. 0 ) THEN 372 idvar = iom_varid( numbdyt, 'vosaline' ) 373 nblendta(igrd) = iom_file(numbdyt)%dimsz(1,idvar) 374 ENDIF 375 WRITE(numout,*) 'Dim size for vosaline is ', nblendta(igrd) 376 ipi = nblendta(igrd) 377 CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 378 379 DO ib = 1, nblen(igrd) 380 DO ik = 1, jpkm1 381 sbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 382 END DO 383 END DO 384 ENDIF ! ln_bdy_tra_frs 385 386 IF(ln_bdy_dyn_frs) THEN 387 388 igrd = 2 ! u-velocity 389 IF ( nblendta(igrd) .le. 0 ) THEN 390 idvar = iom_varid( numbdyu,'vozocrtx' ) 391 nblendta(igrd) = iom_file(numbdyu)%dimsz(1,idvar) 392 ENDIF 393 WRITE(numout,*) 'Dim size for vozocrtx is ', nblendta(igrd) 394 ipi = nblendta(igrd) 395 CALL iom_get ( numbdyu, jpdom_unknown,'vozocrtx',zdta(1:ipi,1:ipj,1:ipk),nbdy_a ) 396 DO ib = 1, nblen(igrd) 397 DO ik = 1, jpkm1 398 ubdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 399 END DO 400 END DO 401 ! 402 igrd = 3 ! v-velocity 403 IF ( nblendta(igrd) .le. 0 ) THEN 404 idvar = iom_varid( numbdyv,'vomecrty' ) 405 nblendta(igrd) = iom_file(numbdyv)%dimsz(1,idvar) 406 ENDIF 407 WRITE(numout,*) 'Dim size for vomecrty is ', nblendta(igrd) 408 ipi = nblendta(igrd) 409 CALL iom_get ( numbdyv, jpdom_unknown,'vomecrty',zdta(1:ipi,1:ipj,1:ipk),nbdy_a ) 410 DO ib = 1, nblen(igrd) 411 DO ik = 1, jpkm1 412 vbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 413 END DO 414 END DO 415 ENDIF ! ln_bdy_dyn_frs 416 417 418 IF ((.NOT.ln_bdy_clim) .AND. (istep(1) > 0)) THEN 419 ! First data time is after start of run 420 ! Put first value in both time levels 421 nbdy_b = nbdy_a 422 IF(ln_bdy_tra_frs) THEN 423 tbdydta(:,:,1) = tbdydta(:,:,2) 424 sbdydta(:,:,1) = sbdydta(:,:,2) 425 ENDIF 426 IF(ln_bdy_dyn_frs) THEN 427 ubdydta(:,:,1) = ubdydta(:,:,2) 428 vbdydta(:,:,1) = vbdydta(:,:,2) 429 ENDIF 430 END IF 431 432 END IF ! nbdy_dta == 0/1 433 434 ! In the case of constant boundary forcing fill bdy arrays once for all 435 IF ((ln_bdy_clim).AND.(ntimes_bdy==1)) THEN 436 IF(ln_bdy_tra_frs) THEN 437 tbdy (:,:) = tbdydta (:,:,2) 438 sbdy (:,:) = sbdydta (:,:,2) 439 ENDIF 440 IF(ln_bdy_dyn_frs) THEN 441 ubdy (:,:) = ubdydta (:,:,2) 442 vbdy (:,:) = vbdydta (:,:,2) 443 ENDIF 444 445 IF(ln_bdy_tra_frs .or. ln_bdy_ice_frs) CALL iom_close( numbdyt ) 446 IF(ln_bdy_dyn_frs) CALL iom_close( numbdyu ) 447 IF(ln_bdy_dyn_frs) CALL iom_close( numbdyv ) 448 END IF 449 450 ENDIF ! End if nit000 451 452 453 ! !---------------------! 454 ! ! at each time step ! 455 ! !---------------------! 456 457 IF( nbdy_dta == 1 .AND. ntimes_bdy > 1 ) THEN 458 ! 459 ! Read one more record if necessary 460 !********************************** 461 462 IF( ln_bdy_clim .AND. imois /= nbdy_b ) THEN ! remember that nbdy_b=0 for kt=nit000 463 nbdy_b = imois 464 nbdy_a = imois + 1 465 nbdy_b = MOD( nbdy_b, iman ) ; IF( nbdy_b == 0 ) nbdy_b = iman 466 nbdy_a = MOD( nbdy_a, iman ) ; IF( nbdy_a == 0 ) nbdy_a = iman 467 lect=.true. 468 ELSEIF( .NOT.ln_bdy_clim .AND. itimer >= istep(nbdy_a) ) THEN 469 470 IF ( nbdy_a < ntimes_bdy ) THEN 471 nbdy_b = nbdy_a 472 nbdy_a = nbdy_a + 1 473 lect =.true. 474 ELSE 475 ! We have reached the end of the file 476 ! put the last data time into both time levels 477 nbdy_b = nbdy_a 478 IF(ln_bdy_tra_frs) THEN 479 tbdydta(:,:,1) = tbdydta(:,:,2) 480 sbdydta(:,:,1) = sbdydta(:,:,2) 481 ENDIF 482 IF(ln_bdy_dyn_frs) THEN 483 ubdydta(:,:,1) = ubdydta(:,:,2) 484 vbdydta(:,:,1) = vbdydta(:,:,2) 485 ENDIF 486 END IF ! nbdy_a < ntimes_bdy 487 488 END IF 489 490 IF( lect ) THEN 491 ! Swap arrays 492 IF(ln_bdy_tra_frs) THEN 493 tbdydta(:,:,1) = tbdydta(:,:,2) 494 sbdydta(:,:,1) = sbdydta(:,:,2) 495 ENDIF 496 IF(ln_bdy_dyn_frs) THEN 497 ubdydta(:,:,1) = ubdydta(:,:,2) 498 vbdydta(:,:,1) = vbdydta(:,:,2) 499 ENDIF 500 501 ! read another set 502 ipj = 1 503 ipk = jpk 504 505 IF(ln_bdy_tra_frs) THEN 506 ! 507 igrd = 1 ! temperature 508 ipi = nblendta(igrd) 509 CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 510 DO ib = 1, nblen(igrd) 511 DO ik = 1, jpkm1 512 tbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 513 END DO 514 END DO 515 ! 516 igrd = 1 ! salinity 517 ipi = nblendta(igrd) 518 CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 519 DO ib = 1, nblen(igrd) 520 DO ik = 1, jpkm1 521 sbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 522 END DO 523 END DO 524 ENDIF ! ln_bdy_tra_frs 525 526 IF(ln_bdy_dyn_frs) THEN 527 ! 528 igrd = 2 ! u-velocity 529 ipi = nblendta(igrd) 530 CALL iom_get ( numbdyu, jpdom_unknown,'vozocrtx',zdta(1:ipi,1:ipj,1:ipk),nbdy_a ) 531 DO ib = 1, nblen(igrd) 532 DO ik = 1, jpkm1 533 ubdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 534 END DO 535 END DO 536 ! 537 igrd = 3 ! v-velocity 538 ipi = nblendta(igrd) 539 CALL iom_get ( numbdyv, jpdom_unknown,'vomecrty',zdta(1:ipi,1:ipj,1:ipk),nbdy_a ) 540 DO ib = 1, nblen(igrd) 541 DO ik = 1, jpkm1 542 vbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 543 END DO 544 END DO 545 ENDIF ! ln_bdy_dyn_frs 546 547 ! 548 IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy_b ',nbdy_b 549 IF(lwp) WRITE(numout,*) '~~~~~~~~ last record file used nbdy_a ',nbdy_a 550 IF (.NOT.ln_bdy_clim) THEN 551 IF(lwp) WRITE(numout,*) 'first record time (s): ', istep(nbdy_b) 552 IF(lwp) WRITE(numout,*) 'model time (s) : ', itimer 553 IF(lwp) WRITE(numout,*) 'second record time (s): ', istep(nbdy_a) 554 ENDIF 555 ! 556 ENDIF ! end lect=.true. 557 558 559 ! Interpolate linearly 560 ! ******************** 561 ! 562 IF( ln_bdy_clim ) THEN ; zxy = FLOAT( nday ) / FLOAT( nobis(nbdy_b) ) + 0.5 - i15 563 ELSE ; zxy = FLOAT( istep(nbdy_b) - itimer ) / FLOAT( istep(nbdy_b) - istep(nbdy_a) ) 564 END IF 565 566 IF(ln_bdy_tra_frs) THEN 567 igrd = 1 ! temperature & salinity 568 DO ib = 1, nblen(igrd) 569 DO ik = 1, jpkm1 570 tbdy(ib,ik) = zxy * tbdydta(ib,ik,2) + (1.-zxy) * tbdydta(ib,ik,1) 571 sbdy(ib,ik) = zxy * sbdydta(ib,ik,2) + (1.-zxy) * sbdydta(ib,ik,1) 572 END DO 573 END DO 574 ENDIF 575 576 IF(ln_bdy_dyn_frs) THEN 577 igrd = 2 ! u-velocity 578 DO ib = 1, nblen(igrd) 579 DO ik = 1, jpkm1 580 ubdy(ib,ik) = zxy * ubdydta(ib,ik,2) + (1.-zxy) * ubdydta(ib,ik,1) 581 END DO 582 END DO 583 ! 584 igrd = 3 ! v-velocity 585 DO ib = 1, nblen(igrd) 586 DO ik = 1, jpkm1 587 vbdy(ib,ik) = zxy * vbdydta(ib,ik,2) + (1.-zxy) * vbdydta(ib,ik,1) 588 END DO 589 END DO 590 ENDIF 591 592 END IF !end if ((nbdy_dta==1).AND.(ntimes_bdy>1)) 593 594 595 ! !---------------------! 596 ! ! last call ! 597 ! !---------------------! 598 IF( kt == nitend ) THEN 599 IF(ln_bdy_tra_frs .or. ln_bdy_ice_frs) CALL iom_close( numbdyt ) ! Closing of the 3 files 600 IF(ln_bdy_dyn_frs) CALL iom_close( numbdyu ) 601 IF(ln_bdy_dyn_frs) CALL iom_close( numbdyv ) 602 ENDIF 603 ! 604 ENDIF ! ln_bdy_dyn_frs .OR. ln_bdy_tra_frs 605 606 END SUBROUTINE bdy_dta 607 608 609 SUBROUTINE bdy_dta_bt( kt, jit ) 69 610 !!--------------------------------------------------------------------------- 70 !! *** SUBROUTINE bdy_dta ***611 !! *** SUBROUTINE bdy_dta_bt *** 71 612 !! 72 !! ** Purpose : Read unstructured boundary data 613 !! ** Purpose : Read unstructured boundary data for Flather condition 73 614 !! 74 !! ** Method : 615 !! ** Method : At the first timestep, read in boundary data for two 616 !! times from the file and time-interpolate. At other 617 !! timesteps, check to see if we need another time from 618 !! the file. If so read it in. Time interpolate. 619 !!--------------------------------------------------------------------------- 620 !!gm DOCTOR names : argument integer : start with "k" 621 INTEGER, INTENT( in ) :: kt ! ocean time-step index 622 INTEGER, INTENT( in ) :: jit ! barotropic time step index 623 ! ! (for timesplitting option, otherwise zero) 75 624 !! 76 !! History : OPA 9.0 ! 05-01 (J. Chanut, A. Sellar) Original77 !! " ! 07-01 (D. Storkey) Update to use IOM module.78 !!---------------------------------------------------------------------------79 !! * Arguments80 INTEGER, INTENT( in ) :: kt ! ocean time-step index81 ! (for timesplitting option, otherwise zero)82 83 !! * Local declarations84 625 LOGICAL :: lect ! flag for reading 85 INTEGER :: jt, jb, jk, jgrd! dummy loop indices626 INTEGER :: it, ib, igrd ! dummy loop indices 86 627 INTEGER :: idvar ! netcdf var ID 87 628 INTEGER :: iman, i15, imois ! Time variables for monthly clim forcing 88 INTEGER :: kbdyt, kbdyu, kbdyv629 INTEGER :: ntimes_bdyt, ntimes_bdyu, ntimes_bdyv 89 630 INTEGER :: itimer, totime 90 631 INTEGER :: ipi, ipj, ipk, inum ! temporary integers (NetCDF read) 91 632 INTEGER :: iyear0, imonth0, iday0 92 633 INTEGER :: ihours0, iminutes0, isec0 93 INTEGER :: kyear, kmonth, kday, ksecs634 INTEGER :: iyear, imonth, iday, isecs 94 635 INTEGER, DIMENSION(jpbtime) :: istept, istepu, istepv ! time arrays from data files 95 REAL(wp) :: dayfrac, zxy, offsett 96 REAL(wp) :: offsetu, offsetv 97 REAL(wp) :: dayjul0, zdayjulini 98 REAL(wp), DIMENSION(jpbtime) :: istepr ! REAL time array from data files 99 REAL(wp), DIMENSION(jpbdta,1,jpk) :: pdta3 ! temporary array for data fields 100 CHARACTER(LEN=80), DIMENSION(3) :: bdyfile 636 REAL(wp) :: dayfrac, zxy, zoffsett 637 REAL(wp) :: zoffsetu, zoffsetv 638 REAL(wp) :: dayjul0, zdayjulini, zntotime 639 REAL(wp) :: zinterval_s, zinterval_e ! First and last interval in time axis 640 REAL(wp), DIMENSION(jpbtime) :: zstepr ! REAL time array from data files 641 REAL(wp), DIMENSION(jpbdta,1) :: zdta ! temporary array for data fields 642 CHARACTER(LEN=80), DIMENSION(3) :: clfile 101 643 CHARACTER(LEN=70 ) :: clunits ! units attribute of time coordinate 102 103 644 !!--------------------------------------------------------------------------- 645 646 !!gm add here the same style as in bdy_dta 647 !!gm clearly bdy_dta_bt and bdy_dta can be combined... 648 !!gm too many things duplicated in the read of data... simplification can be done 104 649 105 650 ! -------------------- ! … … 111 656 ! Some time variables for monthly climatological forcing: 112 657 ! ******************************************************* 658 !!gm here use directely daymod variables 113 659 114 660 iman = INT( raamo ) ! Number of months in a year … … 123 669 ! ********************************************* 124 670 125 itimer = (kt-nit000+1)*rdt ! current time in seconds for interpolation 126 127 ! -------------------- ! 128 ! 1. First call only ! 129 ! -------------------- ! 130 131 IF ( kt == nit000 ) THEN 132 133 istep(:) = 0 134 135 nbdy1=0 136 nbdy2=0 137 138 ! 1.1 Get time information from bdy data file 139 ! ******************************************** 671 itimer = ((kt-1)-nit000+1)*rdt ! current time in seconds for interpolation 672 itimer = itimer + jit*rdtbt ! in non-climatological case 673 674 IF ( ln_bdy_tides ) THEN 675 676 ! -------------------------------------! 677 ! Update BDY fields with tidal forcing ! 678 ! -------------------------------------! 679 680 CALL tide_update( kt, jit ) 681 682 ENDIF 683 684 IF ( ln_bdy_dyn_fla ) THEN 685 686 ! -------------------------------------! 687 ! Update BDY fields with model data ! 688 ! -------------------------------------! 689 690 ! !-------------------! 691 IF( kt == nit000 ) THEN ! First call only ! 692 ! !-------------------! 693 istep_bt(:) = 0 694 nbdy_b_bt = 0 695 nbdy_a_bt = 0 696 697 ! Get time information from bdy data file 698 ! *************************************** 140 699 141 700 IF(lwp) WRITE(numout,*) 142 IF(lwp) WRITE(numout,*) 'bdy_dta :Initialize unstructured boundary data'701 IF(lwp) WRITE(numout,*) 'bdy_dta_bt :Initialize unstructured boundary data for barotropic variables.' 143 702 IF(lwp) WRITE(numout,*) '~~~~~~~' 144 703 145 IF (nbdy_dta == 0) THEN704 IF( nbdy_dta == 0 ) THEN 146 705 IF(lwp) WRITE(numout,*) 'Bdy data are taken from initial conditions' 147 706 … … 154 713 ! necessary time dumps in file are included 155 714 156 bdyfile(1) = filbdy_data_T 157 bdyfile(2) = filbdy_data_U 158 bdyfile(3) = filbdy_data_V 159 160 DO jgrd = 1,3 161 162 CALL iom_open( bdyfile(jgrd), inum ) 163 CALL iom_gettime( inum, istepr, kntime=kbdy, cdunits=clunits ) 164 CALL iom_close( inum ) 715 clfile(1) = filbdy_data_bt_T 716 clfile(2) = filbdy_data_bt_U 717 clfile(3) = filbdy_data_bt_V 718 719 DO igrd = 1,3 720 721 CALL iom_open( clfile(igrd), inum ) 722 CALL iom_gettime( inum, zstepr, kntime=ntimes_bdy, cdunits=clunits ) 723 724 SELECT CASE( igrd ) 725 CASE (1) 726 numbdyt = inum 727 CASE (2) 728 numbdyu = inum 729 CASE (3) 730 numbdyv = inum 731 END SELECT 165 732 166 733 ! Calculate time offset … … 170 737 CALL ymds2ju(iyear0, imonth0, iday0, real(isec0), dayjul0) 171 738 ! Compute model initialization time 172 kyear = ndastp / 10000173 kmonth = ( ndastp - kyear * 10000 ) / 100174 kday = ndastp - kyear * 10000 - kmonth * 100175 ksecs = dayfrac * 86400176 CALL ymds2ju( kyear, kmonth, kday, real(ksecs) , zdayjulini)177 ! offset from initialization date:178 offset = (dayjul0-zdayjulini)*86400739 iyear = ndastp / 10000 740 imonth = ( ndastp - iyear * 10000 ) / 100 741 iday = ndastp - iyear * 10000 - imonth * 100 742 isecs = dayfrac * 86400 743 CALL ymds2ju(iyear, imonth, iday, real(isecs) , zdayjulini) 744 ! zoffset from initialization date: 745 zoffset = (dayjul0-zdayjulini)*86400 179 746 ! 180 747 181 748 7000 FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 182 183 749 184 750 !! TO BE DONE... Check consistency between calendar from file … … 186 752 !! when calendar in model available outside of IOIPSL. 187 753 188 IF(lwp) WRITE(numout,*) 'kdby (number of times): ',kbdy189 IF(lwp) WRITE(numout,*) 'offset: ',offset190 IF(lwp) WRITE(numout,*) 'totime: ',totime191 IF(lwp) WRITE(numout,*) 'istepr: ',istepr192 193 754 ! Check that there are not too many times in the file. 194 IF (kbdy > jpbtime) THEN 195 IF(lwp) WRITE(numout,*) 196 IF(lwp) WRITE(numout,*) ' E R R O R : Number of time dumps in files exceed jpbtime parameter' 197 IF(lwp) WRITE(numout,*) ' ========== jpbtime= ',jpbtime,' kbdy= ', kbdy 198 IF(lwp) WRITE(numout,*) ' Check file: ', bdyfile(jgrd) 199 nstop = nstop + 1 755 IF (ntimes_bdy_bt > jpbtime) CALL ctl_stop( & 756 'Number of time dumps in bdy file exceed jpbtime parameter', & 757 'Check file:' // TRIM(clfile(igrd)) ) 758 759 ! Check that time array increases (or interp will fail): 760 DO it = 2, ntimes_bdy 761 IF ( zstepr(it-1) >= zstepr(it) ) THEN 762 CALL ctl_stop('Time array in unstructured boundary data file', & 763 'does not continuously increase.', & 764 'Check file:' // TRIM(clfile(igrd)) ) 765 EXIT 766 END IF 767 END DO 768 769 IF ( .NOT. ln_bdy_clim ) THEN 770 ! Check that times in file span model run time: 771 772 ! Note: the fields may be time means, so we allow nit000 to be before 773 ! first time in the file, provided that it falls inside the meaning 774 ! period of the first field. Until we can get the meaning period 775 ! from the file, use the interval between fields as a proxy. 776 ! If nit000 is before the first time, use the value at first time 777 ! instead of extrapolating. This is done by putting time 1 into 778 ! both time levels. 779 ! The same applies to the last time level: see setting of lect below. 780 781 IF ( ntimes_bdy == 1 ) CALL ctl_stop( & 782 'There is only one time dump in data files', & 783 'Set ln_bdy_clim=.true. in namelist for constant bdy forcing.' ) 784 785 zinterval_s = zstepr(2) - zstepr(1) 786 zinterval_e = zstepr(ntimes_bdy) - zstepr(ntimes_bdy-1) 787 788 IF ( zstepr(1) - zinterval_s / 2.0 > 0 ) THEN 789 IF(lwp) WRITE(numout,*) 'First bdy time relative to nit000:', zstepr(1) 790 IF(lwp) WRITE(numout,*) 'Interval between first two times: ', zinterval_s 791 CALL ctl_stop( 'First data time is after start of run', & 792 'by more than half a meaning period', & 793 'Check file: ' // TRIM(clfile(igrd)) ) 794 END IF 795 796 IF ( zstepr(ntimes_bdy) + zinterval_e / 2.0 < zntotime ) THEN 797 IF(lwp) WRITE(numout,*) 'Last bdy time relative to nit000:', zstepr(ntimes_bdy) 798 IF(lwp) WRITE(numout,*) 'Interval between last two times: ', zinterval_e 799 CALL ctl_stop( 'Last data time is before end of run', & 800 'by more than half a meaning period', & 801 'Check file: ' // TRIM(clfile(igrd)) ) 802 END IF 803 804 END IF ! .NOT. ln_bdy_clim 805 806 IF ( igrd .EQ. 1) THEN 807 ntimes_bdyt = ntimes_bdy_bt 808 zoffsett = zoffset 809 istept(:) = INT( zstepr(:) + zoffset ) 810 ELSE IF (igrd .EQ. 2) THEN 811 ntimes_bdyu = ntimes_bdy_bt 812 zoffsetu = zoffset 813 istepu(:) = INT( zstepr(:) + zoffset ) 814 ELSE IF (igrd .EQ. 3) THEN 815 ntimes_bdyv = ntimes_bdy_bt 816 zoffsetv = zoffset 817 istepv(:) = INT( zstepr(:) + zoffset ) 200 818 ENDIF 201 819 202 ! Check that time array increases:203 204 jt=1205 DO WHILE ((istepr(jt+1).GT.istepr(jt)).AND.(jt.LE.(kbdy-1)))206 jt=jt+1207 END DO208 209 IF ((jt.NE.kbdy).AND.(kbdy.GT.1)) THEN210 IF(lwp) WRITE(numout,*)211 IF(lwp) WRITE(numout,*) ' E R R O R : Time array in unstructured boundary data files'212 IF(lwp) WRITE(numout,*) ' ========== does not continuously increase. Check file: '213 IF(lwp) WRITE(numout,*) bdyfile(jgrd)214 IF(lwp) WRITE(numout,*)215 nstop = nstop + 1216 ENDIF217 218 ! Check that times in file span model run time:219 220 IF (istepr(1)+offset > 0) THEN221 IF(lwp) WRITE(numout,*)222 IF(lwp) WRITE(numout,*) ' E R R O R : First time dump in bdy file is after model initial time'223 IF(lwp) WRITE(numout,*) ' ========== Check file: ',bdyfile(jgrd)224 IF(lwp) WRITE(numout,*)225 nstop = nstop + 1226 END IF227 228 IF (istepr(kbdy)+offset < totime) THEN229 IF(lwp) WRITE(numout,*)230 IF(lwp) WRITE(numout,*) ' E R R O R : Last time dump in bdy file is before model final time'231 IF(lwp) WRITE(numout,*) ' ========== Check file: ',bdyfile(jgrd)232 IF(lwp) WRITE(numout,*)233 nstop = nstop + 1234 END IF235 236 IF ( jgrd .EQ. 1) THEN237 kbdyt = kbdy238 offsett = offset239 istept(:) = INT( istepr(:) + offset )240 ELSE IF (jgrd .EQ. 2) THEN241 kbdyu = kbdy242 offsetu = offset243 istepu(:) = INT( istepr(:) + offset )244 ELSE IF (jgrd .EQ. 3) THEN245 kbdyv = kbdy246 offsetv = offset247 istepv(:) = INT( istepr(:) + offset )248 ENDIF249 250 820 ENDDO 251 821 252 822 ! Verify time consistency between files 253 823 254 IF (kbdyu.NE.kbdyt .OR. kbdyv.NE.kbdyt) THEN 255 IF(lwp) WRITE(numout,*) 'Bdy data files must have the same number of time dumps' 256 IF(lwp) WRITE(numout,*) 'Multiple time frequencies not implemented yet' 257 nstop = nstop + 1 824 IF ( ntimes_bdyu /= ntimes_bdyt .OR. ntimes_bdyv /= ntimes_bdyt ) THEN 825 CALL ctl_stop( & 826 'Time axis lengths differ between bdy data files', & 827 'Multiple time frequencies not implemented yet' ) 828 ELSE 829 ntimes_bdy_bt = ntimes_bdyt 258 830 ENDIF 259 kbdy = kbdyt 260 261 IF (offsetu.NE.offsett .OR. offsetv.NE.offsett) THEN 262 IF(lwp) WRITE(numout,*) 'Bdy data files must have the same time origin' 263 IF(lwp) WRITE(numout,*) 'Multiple time frequencies not implemented yet' 264 nstop = nstop + 1 831 832 IF (zoffsetu.NE.zoffsett .OR. zoffsetv.NE.zoffsett) THEN 833 CALL ctl_stop( & 834 'Bdy data files must have the same time origin', & 835 'Multiple time frequencies not implemented yet' ) 265 836 ENDIF 266 offset =offsett837 zoffset = zoffsett 267 838 268 839 !! Check that times are the same in the three files... HERE. 269 istep (:) = istept(:)840 istep_bt(:) = istept(:) 270 841 271 842 ! Check number of time dumps: 272 IF ((kbdy.EQ.1).AND.(.NOT.ln_bdy_clim)) THEN273 IF(lwp) WRITE(numout,*)274 IF(lwp) WRITE(numout,*) ' E R R O R : There is only one time dump in data files'275 IF(lwp) WRITE(numout,*) ' ========== Choose ln_bdy_clim=.true. in namelist for constant bdy forcing.'276 IF(lwp) WRITE(numout,*)277 nstop = nstop + 1278 ENDIF279 280 843 IF (ln_bdy_clim) THEN 281 IF ((kbdy.NE.1).AND.(kbdy.NE.12)) THEN 282 IF(lwp) WRITE(numout,*) 283 IF(lwp) WRITE(numout,*) ' E R R O R : For climatological boundary forcing (ln_bdy_clim=.true.),' 284 IF(lwp) WRITE(numout,*) ' ========== bdy data files must contain 1 or 12 time dumps.' 285 IF(lwp) WRITE(numout,*) 286 nstop = nstop + 1 287 ELSEIF (kbdy.EQ.1 ) THEN 844 SELECT CASE ( ntimes_bdy_bt ) 845 CASE( 1 ) 288 846 IF(lwp) WRITE(numout,*) 289 847 IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files' 290 848 IF(lwp) WRITE(numout,*) 291 ELSEIF (kbdy.EQ.12) THEN849 CASE( 12 ) 292 850 IF(lwp) WRITE(numout,*) 293 851 IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files' 294 852 IF(lwp) WRITE(numout,*) 295 ENDIF 853 CASE DEFAULT 854 CALL ctl_stop( & 855 'For climatological boundary forcing (ln_bdy_clim=.true.),',& 856 'bdy data files must contain 1 or 12 time dumps.' ) 857 END SELECT 296 858 ENDIF 297 859 298 860 ! Find index of first record to read (before first model time). 299 861 300 jt=1301 DO WHILE ( ((istep (jt+1)) <= 0 ).AND.(jt.LE.(kbdy-1)))302 jt=jt+1303 END DO 304 nbdy 1 = jt305 306 WRITE(numout,*) 'Time offset is ', offset307 WRITE(numout,*) 'First record to read is ',nbdy 1862 it=1 863 DO WHILE ( ((istep_bt(it+1)) <= 0 ).AND.(it.LE.(ntimes_bdy_bt-1))) 864 it=it+1 865 END DO 866 nbdy_b_bt = it 867 868 WRITE(numout,*) 'Time offset is ',zoffset 869 WRITE(numout,*) 'First record to read is ',nbdy_b_bt 308 870 309 871 ENDIF ! endif (nbdy_dta == 1) … … 314 876 IF ( nbdy_dta == 0) THEN 315 877 ! boundary data arrays are filled with initial conditions 316 DO jk = 1, jpkm1 317 318 #if ! defined key_barotropic 319 jgrd = 1 ! T-points data 320 DO jb = 1, nblen(jgrd) 321 tbdy(jb,jk) = tn(nbi(jb,jgrd), nbj(jb,jgrd), jk) 322 sbdy(jb,jk) = sn(nbi(jb,jgrd), nbj(jb,jgrd), jk) 323 END DO 324 #endif 325 326 jgrd = 2 ! U-points data 327 DO jb = 1, nblen(jgrd) 328 ubdy(jb,jk) = un(nbi(jb,jgrd), nbj(jb,jgrd), jk) 329 END DO 330 331 jgrd = 3 ! V-points data 332 DO jb = 1, nblen(jgrd) 333 vbdy(jb,jk) = vn(nbi(jb,jgrd), nbj(jb,jgrd), jk) 334 END DO 335 END DO 878 igrd = 2 ! U-points data 879 DO ib = 1, nblen(igrd) 880 ubtbdy(ib) = un(nbi(ib,igrd), nbj(ib,igrd), 1) 881 END DO 882 883 igrd = 3 ! V-points data 884 DO ib = 1, nblen(igrd) 885 vbtbdy(ib) = vn(nbi(ib,igrd), nbj(ib,igrd), 1) 886 END DO 887 888 igrd = 1 ! T-points data 889 DO ib = 1, nblen(igrd) 890 sshbdy(ib) = sshn(nbi(ib,igrd), nbj(ib,igrd)) 891 END DO 336 892 337 893 ELSEIF (nbdy_dta == 1) THEN 338 894 339 895 ! Set first record in the climatological case: 340 IF ((ln_bdy_clim).AND.( kbdy==1)) THEN341 nbdy 2= 1342 ELSEIF ((ln_bdy_clim).AND.( kbdy==iman)) THEN343 nbdy 1= 0344 nbdy 2= imois896 IF ((ln_bdy_clim).AND.(ntimes_bdy_bt==1)) THEN 897 nbdy_a_bt = 1 898 ELSEIF ((ln_bdy_clim).AND.(ntimes_bdy_bt==iman)) THEN 899 nbdy_b_bt = 0 900 nbdy_a_bt = imois 345 901 ELSE 346 nbdy 2 = nbdy1902 nbdy_a_bt = nbdy_b_bt 347 903 END IF 348 904 349 905 ! Open Netcdf files: 350 906 351 CALL iom_open ( filbdy_data_ T, bdynumt )352 CALL iom_open ( filbdy_data_ U, bdynumu)353 CALL iom_open ( filbdy_data_ V, bdynumv)907 CALL iom_open ( filbdy_data_bt_T, numbdyt_bt ) 908 CALL iom_open ( filbdy_data_bt_U, numbdyu_bt ) 909 CALL iom_open ( filbdy_data_bt_V, numbdyv_bt ) 354 910 355 911 ! Read first record: 356 912 ipj=1 357 ipk=jpk 358 jgrd=1 359 ipi=nblendta(jgrd) 360 361 #if ! defined key_barotropic 362 !temperature 363 364 jgrd=1 365 IF ( nblendta(jgrd) .le. 0 ) THEN 366 idvar = iom_varid( bdynumt,'votemper' ) 367 nblendta(jgrd) = iom_file(bdynumt)%dimsz(1,idvar) 913 igrd=1 914 ipi=nblendta(igrd) 915 916 ! ssh 917 igrd=1 918 IF ( nblendta(igrd) .le. 0 ) THEN 919 idvar = iom_varid( numbdyt_bt,'sossheig' ) 920 nblendta(igrd) = iom_file(numbdyt_bt)%dimsz(1,idvar) 368 921 ENDIF 369 WRITE(numout,*) 'Dim size for votemper is ',nblendta(jgrd) 370 ipi=nblendta(jgrd) 371 372 CALL iom_get ( bdynumt, jpdom_unknown,'votemper',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 373 374 DO jb=1, nblen(jgrd) 375 DO jk=1,jpkm1 376 tbdydta(jb,jk,2) = pdta3(nbmap(jb,jgrd),1,jk) 377 END DO 378 END DO 379 380 ! salinity 381 382 jgrd=1 383 IF ( nblendta(jgrd) .le. 0 ) THEN 384 idvar = iom_varid( bdynumt,'vosaline' ) 385 nblendta(jgrd) = iom_file(bdynumt)%dimsz(1,idvar) 922 WRITE(numout,*) 'Dim size for sossheig is ',nblendta(igrd) 923 ipi=nblendta(igrd) 924 925 CALL iom_get ( numbdyt_bt, jpdom_unknown,'sossheig',zdta(1:ipi,1:ipj),nbdy_a_bt ) 926 927 DO ib=1, nblen(igrd) 928 sshbdydta(ib,2) = zdta(nbmap(ib,igrd),1) 929 END DO 930 931 ! u-velocity 932 igrd=2 933 IF ( nblendta(igrd) .le. 0 ) THEN 934 idvar = iom_varid( numbdyu_bt,'vobtcrtx' ) 935 nblendta(igrd) = iom_file(numbdyu_bt)%dimsz(1,idvar) 386 936 ENDIF 387 WRITE(numout,*) 'Dim size for vosaline is ',nblendta(jgrd) 388 ipi=nblendta(jgrd) 389 390 CALL iom_get ( bdynumt, jpdom_unknown,'vosaline',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 391 392 DO jb=1, nblen(jgrd) 393 DO jk=1,jpkm1 394 sbdydta(jb,jk,2) = pdta3(nbmap(jb,jgrd),1,jk) 395 END DO 396 END DO 397 #endif 398 399 ! u-velocity 400 401 jgrd=2 402 IF ( nblendta(jgrd) .le. 0 ) THEN 403 idvar = iom_varid( bdynumu,'vozocrtx' ) 404 nblendta(jgrd) = iom_file(bdynumu)%dimsz(1,idvar) 937 WRITE(numout,*) 'Dim size for vobtcrtx is ',nblendta(igrd) 938 ipi=nblendta(igrd) 939 940 CALL iom_get ( numbdyu_bt, jpdom_unknown,'vobtcrtx',zdta(1:ipi,1:ipj),nbdy_a_bt ) 941 942 DO ib=1, nblen(igrd) 943 ubtbdydta(ib,2) = zdta(nbmap(ib,igrd),1) 944 END DO 945 946 ! v-velocity 947 igrd=3 948 IF ( nblendta(igrd) .le. 0 ) THEN 949 idvar = iom_varid( numbdyv_bt,'vobtcrty' ) 950 nblendta(igrd) = iom_file(numbdyv_bt)%dimsz(1,idvar) 405 951 ENDIF 406 WRITE(numout,*) 'Dim size for vozocrtx is ',nblendta(jgrd) 407 ipi=nblendta(jgrd) 408 409 CALL iom_get ( bdynumu, jpdom_unknown,'vozocrtx',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 410 411 DO jb=1, nblen(jgrd) 412 DO jk=1,jpkm1 413 ubdydta(jb,jk,2) = pdta3(nbmap(jb,jgrd),1,jk) 414 END DO 415 END DO 416 417 ! v-velocity 418 jgrd=3 419 IF ( nblendta(jgrd) .le. 0 ) THEN 420 idvar = iom_varid( bdynumv,'vomecrty' ) 421 nblendta(jgrd) = iom_file(bdynumv)%dimsz(1,idvar) 422 ENDIF 423 WRITE(numout,*) 'Dim size for vomecrty is ',nblendta(jgrd) 424 ipi=nblendta(jgrd) 425 426 CALL iom_get ( bdynumv, jpdom_unknown,'vomecrty',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 427 428 DO jb=1, nblen(jgrd) 429 DO jk=1,jpkm1 430 vbdydta(jb,jk,2) = pdta3(nbmap(jb,jgrd),1,jk) 431 END DO 952 WRITE(numout,*) 'Dim size for vobtcrty is ',nblendta(igrd) 953 ipi=nblendta(igrd) 954 955 CALL iom_get ( numbdyv_bt, jpdom_unknown,'vobtcrty',zdta(1:ipi,1:ipj),nbdy_a_bt ) 956 957 DO ib=1, nblen(igrd) 958 vbtbdydta(ib,2) = zdta(nbmap(ib,igrd),1) 432 959 END DO 433 960 … … 435 962 436 963 ! In the case of constant boundary forcing fill bdy arrays once for all 437 IF ((ln_bdy_clim).AND.(kbdy==1)) THEN 438 #if ! defined key_barotropic 439 tbdy (:,:) = tbdydta (:,:,2) 440 sbdy (:,:) = sbdydta (:,:,2) 441 #endif 442 ubdy (:,:) = ubdydta (:,:,2) 443 vbdy (:,:) = vbdydta (:,:,2) 444 445 CALL iom_close( bdynumt ) 446 CALL iom_close( bdynumu ) 447 CALL iom_close( bdynumv ) 964 IF ((ln_bdy_clim).AND.(ntimes_bdy_bt==1)) THEN 965 966 ubtbdy (:) = ubtbdydta (:,2) 967 vbtbdy (:) = vbtbdydta (:,2) 968 sshbdy (:) = sshbdydta (:,2) 969 970 CALL iom_close( numbdyt_bt ) 971 CALL iom_close( numbdyu_bt ) 972 CALL iom_close( numbdyv_bt ) 973 448 974 END IF 449 975 … … 454 980 ! -------------------- ! 455 981 456 IF ((nbdy_dta==1).AND.( kbdy>1)) THEN982 IF ((nbdy_dta==1).AND.(ntimes_bdy_bt>1)) THEN 457 983 458 984 ! 2.1 Read one more record if necessary 459 985 !************************************** 460 986 461 IF ( (ln_bdy_clim).AND.(imois/=nbdy 1) ) THEN ! remember that nbdy1=0 for kt=nit000462 nbdy 1= imois463 nbdy 2= imois+1464 nbdy 1 = MOD( nbdy1, iman )465 IF( nbdy 1 == 0 ) nbdy1= iman466 nbdy 2 = MOD( nbdy2, iman )467 IF( nbdy 2 == 0 ) nbdy2= iman987 IF ( (ln_bdy_clim).AND.(imois/=nbdy_b_bt) ) THEN ! remember that nbdy_b_bt=0 for kt=nit000 988 nbdy_b_bt = imois 989 nbdy_a_bt = imois+1 990 nbdy_b_bt = MOD( nbdy_b_bt, iman ) 991 IF( nbdy_b_bt == 0 ) nbdy_b_bt = iman 992 nbdy_a_bt = MOD( nbdy_a_bt, iman ) 993 IF( nbdy_a_bt == 0 ) nbdy_a_bt = iman 468 994 lect=.true. 469 995 470 ELSEIF ((.NOT.ln_bdy_clim).AND.(itimer >= istep(nbdy2))) THEN 471 nbdy1=nbdy2 472 nbdy2=nbdy2+1 473 lect=.true. 474 END IF 475 476 IF (lect) THEN 477 478 ! Swap arrays 479 #if ! defined key_barotropic 480 tbdydta(:,:,1) = tbdydta(:,:,2) 481 sbdydta(:,:,1) = sbdydta(:,:,2) 482 #endif 483 ubdydta(:,:,1) = ubdydta(:,:,2) 484 vbdydta(:,:,1) = vbdydta(:,:,2) 485 486 ! read another set 487 488 ipj=1 489 ipk=jpk 490 jgrd=1 491 ipi=nblendta(jgrd) 492 493 #if ! defined key_barotropic 494 !temperature 495 496 jgrd=1 497 ipi=nblendta(jgrd) 498 499 CALL iom_get ( bdynumt, jpdom_unknown,'votemper',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 500 501 DO jb=1, nblen(jgrd) 502 DO jk=1,jpkm1 503 tbdydta(jb,jk,2) = pdta3(nbmap(jb,jgrd),1,jk) 504 END DO 505 END DO 506 507 ! salinity 508 509 jgrd=1 510 ipi=nblendta(jgrd) 511 512 CALL iom_get ( bdynumt, jpdom_unknown,'vosaline',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 513 514 DO jb=1, nblen(jgrd) 515 DO jk=1,jpkm1 516 sbdydta(jb,jk,2) = pdta3(nbmap(jb,jgrd),1,jk) 517 END DO 518 END DO 519 #endif 520 521 ! u-velocity 522 523 jgrd=2 524 ipi=nblendta(jgrd) 525 526 CALL iom_get ( bdynumu, jpdom_unknown,'vozocrtx',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 527 528 DO jb=1, nblen(jgrd) 529 DO jk=1,jpkm1 530 ubdydta(jb,jk,2) = pdta3(nbmap(jb,jgrd),1,jk) 531 END DO 532 END DO 533 534 ! v-velocity 535 jgrd=3 536 ipi=nblendta(jgrd) 537 538 CALL iom_get ( bdynumv, jpdom_unknown,'vomecrty',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 539 540 DO jb=1, nblen(jgrd) 541 DO jk=1,jpkm1 542 vbdydta(jb,jk,2) = pdta3(nbmap(jb,jgrd),1,jk) 543 END DO 544 END DO 545 546 IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy1 ',nbdy1 547 IF(lwp) WRITE(numout,*) '~~~~~~~~ last record file used nbdy2 ',nbdy2 548 IF (.NOT.ln_bdy_clim) THEN 549 IF(lwp) WRITE(numout,*) 'first record time (s): ', istep(nbdy1) 550 IF(lwp) WRITE(numout,*) 'model time (s) : ', itimer 551 IF(lwp) WRITE(numout,*) 'second record time (s): ', istep(nbdy2) 552 ENDIF 553 END IF ! end lect=.true. 554 555 556 ! 2.2 Interpolate linearly: 557 ! *************************** 558 559 IF (ln_bdy_clim) THEN 560 zxy = FLOAT( nday ) / FLOAT( nobis(nbdy1) ) + 0.5 - i15 561 ELSE 562 zxy = FLOAT(istep(nbdy1)-itimer) / FLOAT(istep(nbdy1)-istep(nbdy2)) 563 END IF 564 565 #if ! defined key_barotropic 566 jgrd=1 567 DO jb=1, nblen(jgrd) 568 DO jk=1, jpkm1 569 tbdy(jb,jk) = zxy * tbdydta(jb,jk,2) + & 570 (1.-zxy) * tbdydta(jb,jk,1) 571 sbdy(jb,jk) = zxy * sbdydta(jb,jk,2) + & 572 (1.-zxy) * sbdydta(jb,jk,1) 573 END DO 574 END DO 575 #endif 576 577 jgrd=2 578 DO jb=1, nblen(jgrd) 579 DO jk=1, jpkm1 580 ubdy(jb,jk) = zxy * ubdydta(jb,jk,2) + & 581 (1.-zxy) * ubdydta(jb,jk,1) 582 END DO 583 END DO 584 585 jgrd=3 586 DO jb=1, nblen(jgrd) 587 DO jk=1, jpkm1 588 vbdy(jb,jk) = zxy * vbdydta(jb,jk,2) + & 589 (1.-zxy) * vbdydta(jb,jk,1) 590 END DO 591 END DO 592 593 END IF !end if ((nbdy_dta==1).AND.(kbdy>1)) 594 595 ! ------------------- ! 596 ! Last call kt=nitend ! 597 ! ------------------- ! 598 599 ! Closing of the 3 files 600 IF( kt == nitend ) THEN 601 CALL iom_close( bdynumt ) 602 CALL iom_close( bdynumu ) 603 CALL iom_close( bdynumv ) 604 ENDIF 605 606 END SUBROUTINE bdy_dta 607 608 SUBROUTINE bdy_dta_bt( kt, jit ) 609 !!--------------------------------------------------------------------------- 610 !! *** SUBROUTINE bdy_dta_bt *** 611 !! 612 !! ** Purpose : Read unstructured boundary data for barotropic variables 613 !! 614 !! ** Method : 615 !! 616 !! History : OPA 9.0 ! 07-07 (D. Storkey) Original 617 !!--------------------------------------------------------------------------- 618 !! * Arguments 619 INTEGER, INTENT( in ) :: kt ! ocean time-step index 620 INTEGER, INTENT( in ) :: jit ! barotropic time step index 621 ! (for timesplitting option, otherwise zero) 622 623 !! * Local declarations 624 LOGICAL :: lect ! flag for reading 625 INTEGER :: jt, jb, jgrd ! dummy loop indices 626 INTEGER :: idvar ! netcdf var ID 627 INTEGER :: iman, i15, imois ! Time variables for monthly clim forcing 628 INTEGER :: kbdyt, kbdyu, kbdyv 629 INTEGER :: itimer, totime 630 INTEGER :: ipi, ipj, ipk, inum ! temporary integers (NetCDF read) 631 INTEGER :: iyear0, imonth0, iday0 632 INTEGER :: ihours0, iminutes0, isec0 633 INTEGER :: kyear, kmonth, kday, ksecs 634 INTEGER, DIMENSION(jpbtime) :: istept, istepu, istepv ! time arrays from data files 635 REAL(wp) :: dayfrac, zxy, offsett 636 REAL(wp) :: offsetu, offsetv 637 REAL(wp) :: dayjul0, zdayjulini 638 REAL(wp), DIMENSION(jpbtime) :: istepr ! REAL time array from data files 639 REAL(wp), DIMENSION(jpbdta,1) :: pdta2 ! temporary array for data fields 640 REAL(wp), DIMENSION(jpbdta,1,jpk) :: pdta3 ! temporary array for data fields 641 CHARACTER(LEN=80), DIMENSION(3) :: bdyfile 642 CHARACTER(LEN=70 ) :: clunits ! units attribute of time coordinate 643 644 !!--------------------------------------------------------------------------- 645 646 ! -------------------- ! 647 ! Initialization ! 648 ! -------------------- ! 649 650 lect = .false. ! If true, read a time record 651 652 ! Some time variables for monthly climatological forcing: 653 ! ******************************************************* 654 655 iman = INT( raamo ) ! Number of months in a year 656 657 i15 = INT( 2*FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) ) 658 ! i15=0 if the current day is in the first half of the month, else i15=1 659 660 imois = nmonth + i15 - 1 ! imois is the first month record 661 IF( imois == 0 ) imois = iman 662 663 ! Time variable for non-climatological forcing: 664 ! ********************************************* 665 666 itimer = ((kt-1)-nit000+1)*rdt ! current time in seconds for interpolation 667 itimer = itimer + jit*rdtbt ! in non-climatological case 668 669 ! -------------------------------------! 670 ! Update BDY fields with tidal forcing ! 671 ! -------------------------------------! 672 673 IF ( lk_bdy_tides ) CALL tide_update( kt, jit ) 674 675 IF ( lk_bdy ) THEN 676 677 ! -------------------- ! 678 ! 1. First call only ! 679 ! -------------------- ! 680 681 IF ( kt == nit000 ) THEN 682 683 istep_bt(:) = 0 684 685 nbdy1_bt=0 686 nbdy2_bt=0 687 688 ! 1.1 Get time information from bdy data file 689 ! ******************************************** 690 691 IF(lwp) WRITE(numout,*) 692 IF(lwp) WRITE(numout,*) 'bdy_dta_bt :Initialize unstructured boundary data for barotropic variables.' 693 IF(lwp) WRITE(numout,*) '~~~~~~~' 694 695 IF (nbdy_dta == 0) THEN 696 IF(lwp) WRITE(numout,*) 'Bdy data are taken from initial conditions' 697 698 ELSEIF (nbdy_dta == 1) THEN 699 IF(lwp) WRITE(numout,*) 'Bdy data are read in netcdf files' 700 701 dayfrac = adatrj-float(itimer)/86400. ! day fraction at time step kt-1 702 dayfrac = dayfrac - INT(dayfrac) ! 703 totime = (nitend-nit000+1)*rdt ! Total time of the run to verify that all the 704 ! necessary time dumps in file are included 705 706 bdyfile(1) = filbdy_data_bt_T 707 bdyfile(2) = filbdy_data_bt_U 708 bdyfile(3) = filbdy_data_bt_V 709 710 DO jgrd = 1,3 711 712 CALL iom_open( bdyfile(jgrd), inum ) 713 CALL iom_gettime( inum, istepr, kntime=kbdy, cdunits=clunits ) 714 CALL iom_close( inum ) 715 716 ! Calculate time offset 717 READ(clunits,7000) iyear0, imonth0, iday0, ihours0, iminutes0, isec0 718 ! Convert time origin in file to julian days 719 isec0 = isec0 + ihours0*60.*60. + iminutes0*60. 720 CALL ymds2ju(iyear0, imonth0, iday0, real(isec0), dayjul0) 721 ! Compute model initialization time 722 kyear = ndastp / 10000 723 kmonth = ( ndastp - kyear * 10000 ) / 100 724 kday = ndastp - kyear * 10000 - kmonth * 100 725 ksecs = dayfrac * 86400 726 CALL ymds2ju(kyear, kmonth, kday, real(ksecs) , zdayjulini) 727 ! offset from initialization date: 728 offset = (dayjul0-zdayjulini)*86400 729 ! 730 731 7000 FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 732 733 !! TO BE DONE... Check consistency between calendar from file 734 !! (available optionally from iom_gettime) and calendar in model 735 !! when calendar in model available outside of IOIPSL. 736 737 IF(lwp) WRITE(numout,*) 'kdby (number of times): ',kbdy_bt 738 IF(lwp) WRITE(numout,*) 'offset: ',offset 739 IF(lwp) WRITE(numout,*) 'totime: ',totime 740 IF(lwp) WRITE(numout,*) 'istepr: ',istepr 741 742 ! Check that there are not too many times in the file. 743 IF (kbdy_bt > jpbtime) THEN 744 IF(lwp) WRITE(numout,*) 745 IF(lwp) WRITE(numout,*) ' E R R O R : Number of time dumps in files exceed jpbtime parameter' 746 IF(lwp) WRITE(numout,*) ' ========== jpbtime= ',jpbtime,' kbdy_bt= ', kbdy_bt 747 IF(lwp) WRITE(numout,*) ' Check file: ', bdyfile(jgrd) 748 nstop = nstop + 1 749 ENDIF 750 751 ! Check that time array increases: 752 753 jt=1 754 DO WHILE ((istepr(jt+1).GT.istepr(jt)).AND.(jt.LE.(kbdy_bt-1))) 755 jt=jt+1 756 END DO 757 758 IF ((jt.NE.kbdy_bt).AND.(kbdy_bt.GT.1)) THEN 759 IF(lwp) WRITE(numout,*) 760 IF(lwp) WRITE(numout,*) ' E R R O R : Time array in unstructured boundary data files' 761 IF(lwp) WRITE(numout,*) ' ========== does not continuously increase. Check file: ' 762 IF(lwp) WRITE(numout,*) bdyfile(jgrd) 763 IF(lwp) WRITE(numout,*) 764 nstop = nstop + 1 765 ENDIF 766 767 ! Check that times in file span model run time: 768 769 IF (istepr(1)+offset > 0) THEN 770 IF(lwp) WRITE(numout,*) 771 IF(lwp) WRITE(numout,*) ' E R R O R : First time dump in bdy file is after model initial time' 772 IF(lwp) WRITE(numout,*) ' ========== Check file: ',bdyfile(jgrd) 773 IF(lwp) WRITE(numout,*) 774 nstop = nstop + 1 775 END IF 776 777 IF (istepr(kbdy_bt)+offset < totime) THEN 778 IF(lwp) WRITE(numout,*) 779 IF(lwp) WRITE(numout,*) ' E R R O R : Last time dump in bdy file is before model final time' 780 IF(lwp) WRITE(numout,*) ' ========== Check file: ',bdyfile(jgrd) 781 IF(lwp) WRITE(numout,*) 782 nstop = nstop + 1 783 END IF 784 785 IF ( jgrd .EQ. 1) THEN 786 kbdyt = kbdy_bt 787 offsett = offset 788 istept(:) = INT( istepr(:) + offset ) 789 ELSE IF (jgrd .EQ. 2) THEN 790 kbdyu = kbdy_bt 791 offsetu = offset 792 istepu(:) = INT( istepr(:) + offset ) 793 ELSE IF (jgrd .EQ. 3) THEN 794 kbdyv = kbdy_bt 795 offsetv = offset 796 istepv(:) = INT( istepr(:) + offset ) 797 ENDIF 798 799 ENDDO 800 801 ! Verify time consistency between files 802 803 IF (kbdyu.NE.kbdyt .OR. kbdyv.NE.kbdyt) THEN 804 IF(lwp) WRITE(numout,*) 'Bdy data files must have the same number of time dumps' 805 IF(lwp) WRITE(numout,*) 'Multiple time frequencies not implemented yet' 806 nstop = nstop + 1 807 ENDIF 808 kbdy_bt = kbdyt 809 810 IF (offsetu.NE.offsett .OR. offsetv.NE.offsett) THEN 811 IF(lwp) WRITE(numout,*) 'Bdy data files must have the same time origin' 812 IF(lwp) WRITE(numout,*) 'Multiple time frequencies not implemented yet' 813 nstop = nstop + 1 814 ENDIF 815 offset = offsett 816 817 !! Check that times are the same in the three files... HERE. 818 istep_bt(:) = istept(:) 819 820 ! Check number of time dumps: 821 IF ((kbdy_bt.EQ.1).AND.(.NOT.ln_bdy_clim)) THEN 822 IF(lwp) WRITE(numout,*) 823 IF(lwp) WRITE(numout,*) ' E R R O R : There is only one time dump in data files' 824 IF(lwp) WRITE(numout,*) ' ========== Choose ln_bdy_clim=.true. in namelist for constant bdy forcing.' 825 IF(lwp) WRITE(numout,*) 826 nstop = nstop + 1 827 ENDIF 828 829 IF (ln_bdy_clim) THEN 830 IF ((kbdy_bt.NE.1).AND.(kbdy_bt.NE.12)) THEN 831 IF(lwp) WRITE(numout,*) 832 IF(lwp) WRITE(numout,*) ' E R R O R : For climatological boundary forcing (ln_bdy_clim=.true.),' 833 IF(lwp) WRITE(numout,*) ' ========== bdy data files must contain 1 or 12 time dumps.' 834 IF(lwp) WRITE(numout,*) 835 nstop = nstop + 1 836 ELSEIF (kbdy_bt.EQ.1 ) THEN 837 IF(lwp) WRITE(numout,*) 838 IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files' 839 IF(lwp) WRITE(numout,*) 840 ELSEIF (kbdy_bt.EQ.12) THEN 841 IF(lwp) WRITE(numout,*) 842 IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files' 843 IF(lwp) WRITE(numout,*) 844 ENDIF 845 ENDIF 846 847 ! Find index of first record to read (before first model time). 848 849 jt=1 850 DO WHILE ( ((istep_bt(jt+1)) <= 0 ).AND.(jt.LE.(kbdy_bt-1))) 851 jt=jt+1 852 END DO 853 nbdy1_bt = jt 854 855 WRITE(numout,*) 'Time offset is ',offset 856 WRITE(numout,*) 'First record to read is ',nbdy1_bt 857 858 ENDIF ! endif (nbdy_dta == 1) 859 860 ! 1.2 Read first record in file if necessary (ie if nbdy_dta == 1) 861 ! ***************************************************************** 862 863 IF ( nbdy_dta == 0) THEN 864 ! boundary data arrays are filled with initial conditions 865 jgrd = 2 ! U-points data 866 DO jb = 1, nblen(jgrd) 867 ubtbdy(jb) = un(nbi(jb,jgrd), nbj(jb,jgrd), 1) 868 END DO 869 870 jgrd = 3 ! V-points data 871 DO jb = 1, nblen(jgrd) 872 vbtbdy(jb) = vn(nbi(jb,jgrd), nbj(jb,jgrd), 1) 873 END DO 874 875 jgrd = 1 ! T-points data 876 DO jb = 1, nblen(jgrd) 877 sshbdy(jb) = sshn(nbi(jb,jgrd), nbj(jb,jgrd)) 878 END DO 879 880 ELSEIF (nbdy_dta == 1) THEN 881 882 ! Set first record in the climatological case: 883 IF ((ln_bdy_clim).AND.(kbdy_bt==1)) THEN 884 nbdy2_bt = 1 885 ELSEIF ((ln_bdy_clim).AND.(kbdy_bt==iman)) THEN 886 nbdy1_bt = 0 887 nbdy2_bt = imois 888 ELSE 889 nbdy2_bt = nbdy1_bt 890 END IF 891 892 ! Open Netcdf files: 893 894 CALL iom_open ( filbdy_data_bt_T, bdynumt_bt ) 895 CALL iom_open ( filbdy_data_bt_U, bdynumu_bt ) 896 CALL iom_open ( filbdy_data_bt_V, bdynumv_bt ) 897 898 ! Read first record: 899 ipj=1 900 jgrd=1 901 ipi=nblendta(jgrd) 902 903 ! ssh 904 jgrd=1 905 IF ( nblendta(jgrd) .le. 0 ) THEN 906 idvar = iom_varid( bdynumt_bt,'sossheig' ) 907 nblendta(jgrd) = iom_file(bdynumt_bt)%dimsz(1,idvar) 908 ENDIF 909 WRITE(numout,*) 'Dim size for sossheig is ',nblendta(jgrd) 910 ipi=nblendta(jgrd) 911 912 CALL iom_get ( bdynumt_bt, jpdom_unknown,'sossheig',pdta2(1:ipi,1:ipj),nbdy2_bt ) 913 914 DO jb=1, nblen(jgrd) 915 sshbdydta(jb,2) = pdta2(nbmap(jb,jgrd),1) 916 END DO 917 918 ! u-velocity 919 jgrd=2 920 IF ( nblendta(jgrd) .le. 0 ) THEN 921 idvar = iom_varid( bdynumu_bt,'vobtcrtx' ) 922 nblendta(jgrd) = iom_file(bdynumu_bt)%dimsz(1,idvar) 923 ENDIF 924 WRITE(numout,*) 'Dim size for vobtcrtx is ',nblendta(jgrd) 925 ipi=nblendta(jgrd) 926 927 CALL iom_get ( bdynumu_bt, jpdom_unknown,'vobtcrtx',pdta2(1:ipi,1:ipj),nbdy2_bt ) 928 929 DO jb=1, nblen(jgrd) 930 ubtbdydta(jb,2) = pdta2(nbmap(jb,jgrd),1) 931 END DO 932 933 ! v-velocity 934 jgrd=3 935 IF ( nblendta(jgrd) .le. 0 ) THEN 936 idvar = iom_varid( bdynumv_bt,'vobtcrty' ) 937 nblendta(jgrd) = iom_file(bdynumv_bt)%dimsz(1,idvar) 938 ENDIF 939 WRITE(numout,*) 'Dim size for vobtcrty is ',nblendta(jgrd) 940 ipi=nblendta(jgrd) 941 942 CALL iom_get ( bdynumv_bt, jpdom_unknown,'vobtcrty',pdta2(1:ipi,1:ipj),nbdy2_bt ) 943 944 DO jb=1, nblen(jgrd) 945 vbtbdydta(jb,2) = pdta2(nbmap(jb,jgrd),1) 946 END DO 947 948 END IF 949 950 ! In the case of constant boundary forcing fill bdy arrays once for all 951 IF ((ln_bdy_clim).AND.(kbdy_bt==1)) THEN 952 953 ubtbdy (:) = ubtbdydta (:,2) 954 vbtbdy (:) = vbtbdydta (:,2) 955 sshbdy (:) = sshbdydta (:,2) 956 957 CALL iom_close( bdynumt_bt ) 958 CALL iom_close( bdynumu_bt ) 959 CALL iom_close( bdynumv_bt ) 960 961 END IF 962 963 ENDIF ! End if nit000 964 965 ! -------------------- ! 966 ! 2. At each time step ! 967 ! -------------------- ! 968 969 IF ((nbdy_dta==1).AND.(kbdy_bt>1)) THEN 970 971 ! 2.1 Read one more record if necessary 972 !************************************** 973 974 IF ( (ln_bdy_clim).AND.(imois/=nbdy1_bt) ) THEN ! remember that nbdy1_bt=0 for kt=nit000 975 nbdy1_bt = imois 976 nbdy2_bt = imois+1 977 nbdy1_bt = MOD( nbdy1_bt, iman ) 978 IF( nbdy1_bt == 0 ) nbdy1_bt = iman 979 nbdy2_bt = MOD( nbdy2_bt, iman ) 980 IF( nbdy2_bt == 0 ) nbdy2_bt = iman 981 lect=.true. 982 983 ELSEIF ((.NOT.ln_bdy_clim).AND.(itimer >= istep_bt(nbdy2_bt))) THEN 984 nbdy1_bt=nbdy2_bt 985 nbdy2_bt=nbdy2_bt+1 996 ELSEIF ((.NOT.ln_bdy_clim).AND.(itimer >= istep_bt(nbdy_a_bt))) THEN 997 nbdy_b_bt=nbdy_a_bt 998 nbdy_a_bt=nbdy_a_bt+1 986 999 lect=.true. 987 1000 END IF … … 998 1011 ipj=1 999 1012 ipk=jpk 1000 jgrd=11001 ipi=nblendta( jgrd)1013 igrd=1 1014 ipi=nblendta(igrd) 1002 1015 1003 1016 1004 1017 ! ssh 1005 jgrd=11006 ipi=nblendta( jgrd)1007 1008 CALL iom_get ( bdynumt_bt, jpdom_unknown,'sossheig',pdta2(1:ipi,1:ipj),nbdy2_bt )1009 1010 DO jb=1, nblen(jgrd)1011 sshbdydta( jb,2) = pdta2(nbmap(jb,jgrd),1)1018 igrd=1 1019 ipi=nblendta(igrd) 1020 1021 CALL iom_get ( numbdyt_bt, jpdom_unknown,'sossheig',zdta(1:ipi,1:ipj),nbdy_a_bt ) 1022 1023 DO ib=1, nblen(igrd) 1024 sshbdydta(ib,2) = zdta(nbmap(ib,igrd),1) 1012 1025 END DO 1013 1026 1014 1027 ! u-velocity 1015 jgrd=21016 ipi=nblendta( jgrd)1017 1018 CALL iom_get ( bdynumu_bt, jpdom_unknown,'vobtcrtx',pdta2(1:ipi,1:ipj),nbdy2_bt )1019 1020 DO jb=1, nblen(jgrd)1021 ubtbdydta( jb,2) = pdta3(nbmap(jb,jgrd),1)1028 igrd=2 1029 ipi=nblendta(igrd) 1030 1031 CALL iom_get ( numbdyu_bt, jpdom_unknown,'vobtcrtx',zdta(1:ipi,1:ipj),nbdy_a_bt ) 1032 1033 DO ib=1, nblen(igrd) 1034 ubtbdydta(ib,2) = zdta(nbmap(ib,igrd),1) 1022 1035 END DO 1023 1036 1024 1037 ! v-velocity 1025 jgrd=31026 ipi=nblendta( jgrd)1027 1028 CALL iom_get ( bdynumv_bt, jpdom_unknown,'vobtcrty',pdta2(1:ipi,1:ipj),nbdy2_bt )1029 1030 DO jb=1, nblen(jgrd)1031 vbtbdydta( jb,2) = pdta3(nbmap(jb,jgrd),1)1032 END DO 1033 1034 1035 IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy 1_bt ',nbdy1_bt1036 IF(lwp) WRITE(numout,*) '~~~~~~~~ last record file used nbdy 2_bt ',nbdy2_bt1038 igrd=3 1039 ipi=nblendta(igrd) 1040 1041 CALL iom_get ( numbdyv_bt, jpdom_unknown,'vobtcrty',zdta(1:ipi,1:ipj),nbdy_a_bt ) 1042 1043 DO ib=1, nblen(igrd) 1044 vbtbdydta(ib,2) = zdta(nbmap(ib,igrd),1) 1045 END DO 1046 1047 1048 IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy_b_bt ',nbdy_b_bt 1049 IF(lwp) WRITE(numout,*) '~~~~~~~~ last record file used nbdy_a_bt ',nbdy_a_bt 1037 1050 IF (.NOT.ln_bdy_clim) THEN 1038 IF(lwp) WRITE(numout,*) 'first record time (s): ', istep_bt(nbdy 1_bt)1051 IF(lwp) WRITE(numout,*) 'first record time (s): ', istep_bt(nbdy_b_bt) 1039 1052 IF(lwp) WRITE(numout,*) 'model time (s) : ', itimer 1040 IF(lwp) WRITE(numout,*) 'second record time (s): ', istep_bt(nbdy 2_bt)1053 IF(lwp) WRITE(numout,*) 'second record time (s): ', istep_bt(nbdy_a_bt) 1041 1054 ENDIF 1042 1055 END IF ! end lect=.true. … … 1047 1060 1048 1061 IF (ln_bdy_clim) THEN 1049 zxy = FLOAT( nday ) / FLOAT( nobis(nbdy 1_bt) ) + 0.5 - i151062 zxy = FLOAT( nday ) / FLOAT( nobis(nbdy_b_bt) ) + 0.5 - i15 1050 1063 ELSE 1051 zxy = FLOAT(istep_bt(nbdy 1_bt)-itimer) / FLOAT(istep_bt(nbdy1_bt)-istep_bt(nbdy2_bt))1064 zxy = FLOAT(istep_bt(nbdy_b_bt)-itimer) / FLOAT(istep_bt(nbdy_b_bt)-istep_bt(nbdy_a_bt)) 1052 1065 END IF 1053 1066 1054 jgrd=11055 DO jb=1, nblen(jgrd)1056 sshbdy( jb) = zxy * sshbdydta(jb,2) + &1057 (1.-zxy) * sshbdydta( jb,1)1058 END DO 1059 1060 jgrd=21061 DO jb=1, nblen(jgrd)1062 ubtbdy( jb) = zxy * ubtbdydta(jb,2) + &1063 (1.-zxy) * ubtbdydta( jb,1)1064 END DO 1065 1066 jgrd=31067 DO jb=1, nblen(jgrd)1068 vbtbdy( jb) = zxy * vbtbdydta(jb,2) + &1069 (1.-zxy) * vbtbdydta( jb,1)1070 END DO 1071 1072 1073 END IF !end if ((nbdy_dta==1).AND.( kbdy_bt>1))1067 igrd=1 1068 DO ib=1, nblen(igrd) 1069 sshbdy(ib) = zxy * sshbdydta(ib,2) + & 1070 (1.-zxy) * sshbdydta(ib,1) 1071 END DO 1072 1073 igrd=2 1074 DO ib=1, nblen(igrd) 1075 ubtbdy(ib) = zxy * ubtbdydta(ib,2) + & 1076 (1.-zxy) * ubtbdydta(ib,1) 1077 END DO 1078 1079 igrd=3 1080 DO ib=1, nblen(igrd) 1081 vbtbdy(ib) = zxy * vbtbdydta(ib,2) + & 1082 (1.-zxy) * vbtbdydta(ib,1) 1083 END DO 1084 1085 1086 END IF !end if ((nbdy_dta==1).AND.(ntimes_bdy_bt>1)) 1074 1087 1075 1088 ! ------------------- ! … … 1079 1092 ! Closing of the 3 files 1080 1093 IF( kt == nitend ) THEN 1081 CALL iom_close( bdynumt_bt )1082 CALL iom_close( bdynumu_bt )1083 CALL iom_close( bdynumv_bt )1094 CALL iom_close( numbdyt_bt ) 1095 CALL iom_close( numbdyu_bt ) 1096 CALL iom_close( numbdyv_bt ) 1084 1097 ENDIF 1085 1098 1086 ENDIF ! l k_bdy1099 ENDIF ! ln_bdy_dyn_frs 1087 1100 1088 1101 END SUBROUTINE bdy_dta_bt … … 1090 1103 1091 1104 #else 1092 !!---------------------------------------------------------------------- --------1093 !! default option: Dummy moduleNO Unstruct Open Boundary Conditions1094 !!---------------------------------------------------------------------- --------1105 !!---------------------------------------------------------------------- 1106 !! Dummy module NO Unstruct Open Boundary Conditions 1107 !!---------------------------------------------------------------------- 1095 1108 CONTAINS 1096 SUBROUTINE bdy_dta( kt ) ! Dummy routine 1097 INTEGER, INTENT (in) :: kt 1109 SUBROUTINE bdy_dta( kt ) ! Empty routine 1098 1110 WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt 1099 1111 END SUBROUTINE bdy_dta 1100 1101 SUBROUTINE bdy_dta_bt( kt, jit ) ! Dummy routine 1102 INTEGER, INTENT (in) :: kt, jit 1103 WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt, jit 1112 SUBROUTINE bdy_dta_bt( kt, kit ) ! Empty routine 1113 WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt, kit 1104 1114 END SUBROUTINE bdy_dta_bt 1105 1115 #endif -
trunk/NEMO/OPA_SRC/BDY/bdydyn.F90
r911 r1125 1 1 MODULE bdydyn 2 !!====================================================================== ===========2 !!====================================================================== 3 3 !! *** MODULE bdydyn *** 4 !! Ocean dynamics: Flow relaxation scheme of velocities on unstruc. open boundary 5 !!================================================================================= 6 7 !!--------------------------------------------------------------------------------- 4 !! Unstructured Open Boundary Cond. : Flow relaxation scheme on velocities 5 !!====================================================================== 6 !! History : 1.0 ! 2005-02 (J. Chanut, A. Sellar) Original code 7 !! - ! 2007-07 (D. Storkey) Move Flather implementation to separate routine. 8 !! 3.0 ! 2008-04 (NEMO team) add in the reference version 9 !!---------------------------------------------------------------------- 10 #if defined key_bdy 11 !!---------------------------------------------------------------------- 12 !! 'key_bdy' : Unstructured Open Boundary Condition 13 !!---------------------------------------------------------------------- 8 14 !! bdy_dyn_frs : relaxation of velocities on unstructured open boundary 9 15 !! bdy_dyn_fla : Flather condition for barotropic solution 10 !!--------------------------------------------------------------------------------- 11 #if defined key_bdy || defined key_bdy_tides 12 !!--------------------------------------------------------------------------------- 13 !! * Modules used 16 !!---------------------------------------------------------------------- 14 17 USE oce ! ocean dynamics and tracers 15 18 USE dom_oce ! ocean space and time domain … … 19 22 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 20 23 USE bdytides ! for tidal harmonic forcing at boundary 21 USE in_out_manager 24 USE in_out_manager ! 22 25 23 26 IMPLICIT NONE 24 27 PRIVATE 25 28 26 !! * Accessibility 27 PUBLIC bdy_dyn_frs ! routine called in dynspg_flt (free surface case ONLY) 28 #if defined key_dynspg_exp || defined key_dynspg_ts 29 PUBLIC bdy_dyn_fla ! routine called in dynspg_flt (free surface case ONLY) 30 #endif 29 PUBLIC bdy_dyn_frs ! routine called in dynspg_flt (free surface case ONLY) 30 # if defined key_dynspg_exp || defined key_dynspg_ts 31 PUBLIC bdy_dyn_fla ! routine called in dynspg_flt (free surface case ONLY) 32 # endif 31 33 32 !!--------------------------------------------------------------------------------- 34 !!---------------------------------------------------------------------- 35 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 36 !! $Id: $ 37 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 38 !!---------------------------------------------------------------------- 33 39 34 40 CONTAINS 35 41 36 SUBROUTINE bdy_dyn_frs 37 !!---------------------------------------------------------------------- --------38 !! SUBROUTINE bdy_dyn_frs39 !! ***********************42 SUBROUTINE bdy_dyn_frs( kt ) 43 !!---------------------------------------------------------------------- 44 !! *** SUBROUTINE bdy_dyn_frs *** 45 !! 40 46 !! ** Purpose : - Apply the Flow Relaxation Scheme for dynamic in the 41 47 !! case of unstructured open boundaries. … … 44 50 !! a three-dimensional baroclinic ocean model with realistic 45 51 !! topography. Tellus, 365-382. 52 !!---------------------------------------------------------------------- 53 INTEGER, INTENT( in ) :: kt ! Main time step counter 46 54 !! 47 !! History :48 !! 9.0 ! 05-02 (J. Chanut, A. Sellar) Original49 !! ! 07-07 (D. Storkey) Move Flather implementation to separate routine.50 !!---------------------------------------------------------------------- --------51 ! ! * Arguments52 I NTEGER, INTENT( in ) :: kt ! Main time step counter55 INTEGER :: ib, ik, igrd ! dummy loop indices 56 INTEGER :: ii, ij ! 2D addresses 57 REAL(wp) :: zwgt ! boundary weight 58 !!---------------------------------------------------------------------- 59 ! 60 IF(ln_bdy_dyn_frs) THEN ! If this is false, then this routine does nothing. 53 61 54 !! * Local declarations 55 INTEGER :: jb, jk, jgrd ! dummy loop indices 56 INTEGER :: ii, ij ! 2D addresses 57 REAL(wp) :: zwgt ! boundary weight 58 59 !!------------------------------------------------------------------------------ 60 !! OPA 9.0, LODYC-IPSL (2003) 61 !!------------------------------------------------------------------------------ 62 63 ! ---------------------------! 64 ! FRS on the total velocity :! 65 ! ---------------------------! 66 67 jgrd=2 !: Relaxation of zonal velocity 68 DO jb = 1, nblen(jgrd) 69 DO jk = 1, jpkm1 70 ii = nbi(jb,jgrd) 71 ij = nbj(jb,jgrd) 72 zwgt = nbw(jb,jgrd) 73 74 ua(ii,ij,jk) = ( ua(ii,ij,jk)*(1.-zwgt) + ubdy(jb,jk)*zwgt ) & 75 * umask(ii,ij,jk) 76 END DO 77 END DO 78 79 jgrd=3 !: Relaxation of meridional velocity 80 DO jb = 1, nblen(jgrd) 81 DO jk = 1, jpkm1 82 ii = nbi(jb,jgrd) 83 ij = nbj(jb,jgrd) 84 zwgt = nbw(jb,jgrd) 85 86 va(ii,ij,jk) = ( va(ii,ij,jk)*(1.-zwgt) + vbdy(jb,jk)*zwgt ) & 87 * vmask(ii,ij,jk) 88 END DO 89 END DO 90 91 CALL lbc_lnk( ua, 'U', 1. ) ! Boundary points should be updated 92 CALL lbc_lnk( va, 'V', 1. ) ! 62 IF( kt == nit000 ) THEN 63 IF(lwp) WRITE(numout,*) 64 IF(lwp) WRITE(numout,*) 'bdy_dyn : Flow Relaxation Scheme on momentum' 65 IF(lwp) WRITE(numout,*) '~~~~~~~' 66 ENDIF 67 ! 68 igrd = 2 ! Relaxation of zonal velocity 69 DO ib = 1, nblen(igrd) 70 DO ik = 1, jpkm1 71 ii = nbi(ib,igrd) 72 ij = nbj(ib,igrd) 73 zwgt = nbw(ib,igrd) 74 ua(ii,ij,ik) = ( ua(ii,ij,ik) * ( 1.- zwgt ) + ubdy(ib,ik) * zwgt ) * umask(ii,ij,ik) 75 END DO 76 END DO 77 ! 78 igrd = 3 ! Relaxation of meridional velocity 79 DO ib = 1, nblen(igrd) 80 DO ik = 1, jpkm1 81 ii = nbi(ib,igrd) 82 ij = nbj(ib,igrd) 83 zwgt = nbw(ib,igrd) 84 va(ii,ij,ik) = ( va(ii,ij,ik) * ( 1.- zwgt ) + vbdy(ib,ik) * zwgt ) * vmask(ii,ij,ik) 85 END DO 86 END DO 87 ! 88 CALL lbc_lnk( ua, 'U', 1. ) ! Boundary points should be updated 89 CALL lbc_lnk( va, 'V', 1. ) ! 90 ! 91 ENDIF ! ln_bdy_dyn_frs 93 92 94 93 END SUBROUTINE bdy_dyn_frs 94 95 95 96 96 #if defined key_dynspg_exp || defined key_dynspg_ts 97 97 !! Option to use Flather with dynspg_flt not coded yet... 98 98 SUBROUTINE bdy_dyn_fla 99 !!---------------------------------------------------------------------- --------100 !! SUBROUTINE bdy_dyn_fla101 !! ***********************99 !!---------------------------------------------------------------------- 100 !! *** SUBROUTINE bdy_dyn_fla *** 101 !! 102 102 !! - Apply Flather boundary conditions on normal barotropic velocities 103 !! (ln_bdy_ fla=.true.)103 !! (ln_bdy_dyn_fla=.true. or ln_bdy_tides=.true.) 104 104 !! 105 105 !! ** WARNINGS about FLATHER implementation: … … 113 113 !! ssh in the code is not updated). 114 114 !! 115 !! - Flather, R. A., 1976: A tidal model of the northwest European 116 !! continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164. 117 !! History : 118 !! 9.0 ! 05-02 (J. Chanut, A. Sellar) Original 119 !! ! 07-07 (D. Storkey) Flather algorithm in separate routine. 120 !!------------------------------------------------------------------------------ 121 !! * Local declarations 122 INTEGER :: jb, jk, jgrd, ji, jj, jim1, jip1, jjm1, jjp1 ! dummy loop indices 123 INTEGER :: ii, ij ! 2D addresses 124 REAL(wp) :: corr ! Flather correction 125 REAL(wp) :: zwgt ! boundary weight 126 REAL(wp) :: elapsed 127 128 !!------------------------------------------------------------------------------ 129 !! OPA 9.0, LODYC-IPSL (2003) 130 !!------------------------------------------------------------------------------ 115 !! References: Flather, R. A., 1976: A tidal model of the northwest European 116 !! continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164. 117 !!---------------------------------------------------------------------- 118 INTEGER :: ib, igrd ! dummy loop indices 119 INTEGER :: ii, ij, iim1, iip1, ijm1, ijp1 ! 2D addresses 120 REAL(wp) :: zcorr ! Flather correction 121 !!---------------------------------------------------------------------- 131 122 132 123 ! ---------------------------------! … … 134 125 ! ---------------------------------! 135 126 127 IF(ln_bdy_dyn_fla .OR. ln_bdy_tides) THEN ! If these are both false, then this routine does nothing. 128 136 129 ! Fill temporary array with ssh data (here spgu): 137 jgrd = 1 138 DO jb = 1, nblenrim(jgrd) 139 ji = nbi(jb,jgrd) 140 jj = nbj(jb,jgrd) 141 spgu(ji, jj) = sshbdy(jb) 142 #if defined key_bdy_tides 143 spgu(ji, jj) = spgu(ji, jj) + sshtide(jb) 144 #endif 130 igrd = 1 131 spgu(:,:) = 0.0 132 DO ib = 1, nblenrim(igrd) 133 ii = nbi(ib,igrd) 134 ij = nbj(ib,igrd) 135 IF( ln_bdy_dyn_fla ) spgu(ii, ij) = sshbdy(ib) 136 IF( ln_bdy_tides ) spgu(ii, ij) = spgu(ii, ij) + sshtide(ib) 145 137 END DO 146 147 jgrd = 2 !: Flather bc on u-velocity; 148 ! remember that flagu=-1 if normal velocity direction is outward 149 ! I think we should rather use after ssh ? 150 151 DO jb = 1, nblenrim(jgrd) 152 ji = nbi(jb,jgrd) 153 jj = nbj(jb,jgrd) 154 jim1 = ji+MAX(0, INT(flagu(jb))) ! T pts i-indice inside the boundary 155 jip1 = ji-MIN(0, INT(flagu(jb))) ! T pts i-indice outside the boundary 156 157 corr = - flagu(jb) * sqrt (grav / (hu_e(ji, jj) + 1.e-20) ) & 158 * ( sshn_e(jim1, jj) - spgu(jip1,jj) ) 159 ua_e(ji, jj) = ( ubtbdy(jb) + utide(jb) ) * hu_e(ji,jj) 160 if ( ln_bdy_fla ) then 161 ua_e(ji,jj) = ua_e(ji,jj) + corr * umask(ji,jj,1) * hu_e(ji,jj) 162 endif 163 138 ! 139 igrd = 2 ! Flather bc on u-velocity; 140 ! ! remember that flagu=-1 if normal velocity direction is outward 141 ! ! I think we should rather use after ssh ? 142 DO ib = 1, nblenrim(igrd) 143 ii = nbi(ib,igrd) 144 ij = nbj(ib,igrd) 145 iim1 = ii + MAX( 0, INT( flagu(ib) ) ) ! T pts i-indice inside the boundary 146 iip1 = ii - MIN( 0, INT( flagu(ib) ) ) ! T pts i-indice outside the boundary 147 ! 148 zcorr = - flagu(ib) * SQRT( grav / (hu_e(ii, ij) + 1.e-20) ) * ( sshn_e(iim1, ij) - spgu(iip1,ij) ) 149 ua_e(ii, ij) = ( ubtbdy(ib) + utide(ib) ) * hu_e(ii,ij) 150 ua_e(ii,ij) = ua_e(ii,ij) + zcorr * umask(ii,ij,1) * hu_e(ii,ij) 164 151 END DO 165 166 jgrd = 3 !: Flather bc on v-velocity 167 ! remember that flagv=-1 if normal velocity direction is outward 168 169 DO jb = 1, nblenrim(jgrd) 170 ji = nbi(jb,jgrd) 171 jj = nbj(jb,jgrd) 172 jjm1 = jj+MAX(0, INT(flagv(jb))) ! T pts j-indice inside the boundary 173 jjp1 = jj-MIN(0, INT(flagv(jb))) ! T pts j-indice outside the boundary 174 175 corr = - flagv(jb) * sqrt (grav / (hv_e(ji, jj) + 1.e-20) ) & 176 * ( sshn_e(ji, jjm1) - spgu(ji,jjp1) ) 177 va_e(ji, jj) = ( vbtbdy(jb) + vtide(jb) ) * hv_e(ji,jj) 178 if ( ln_bdy_fla ) then 179 va_e(ji,jj) = va_e(ji,jj) + corr * vmask(ji,jj,1) * hv_e(ji,jj) 180 endif 181 152 ! 153 igrd = 3 ! Flather bc on v-velocity 154 ! ! remember that flagv=-1 if normal velocity direction is outward 155 DO ib = 1, nblenrim(igrd) 156 ii = nbi(ib,igrd) 157 ij = nbj(ib,igrd) 158 ijm1 = ij + MAX( 0, INT( flagv(ib) ) ) ! T pts j-indice inside the boundary 159 ijp1 = ij - MIN( 0, INT( flagv(ib) ) ) ! T pts j-indice outside the boundary 160 ! 161 zcorr = - flagv(ib) * SQRT( grav / (hv_e(ii, ij) + 1.e-20) ) * ( sshn_e(ii, ijm1) - spgu(ii,ijp1) ) 162 va_e(ii, ij) = ( vbtbdy(ib) + vtide(ib) ) * hv_e(ii,ij) 163 va_e(ii,ij) = va_e(ii,ij) + zcorr * vmask(ii,ij,1) * hv_e(ii,ij) 182 164 END DO 183 165 ! 184 166 CALL lbc_lnk( ua_e, 'U', 1. ) ! Boundary points should be updated 185 167 CALL lbc_lnk( va_e, 'V', 1. ) ! 168 ! 169 ENDIF ! ln_bdy_dyn_fla .or. ln_bdy_tides 186 170 187 171 END SUBROUTINE bdy_dyn_fla … … 189 173 190 174 #else 191 !!================================================================================= 192 !! *** MODULE bdydyn *** 193 !! Ocean dynamics: Flow relaxation scheme of velocities on unstruc. open boundary 194 !!================================================================================= 175 !!---------------------------------------------------------------------- 176 !! Dummy module NO Unstruct Open Boundary Conditions 177 !!---------------------------------------------------------------------- 195 178 CONTAINS 196 197 SUBROUTINE bdy_dyn_frs 198 ! No Unstructured open boundaries ==> empty routine 179 SUBROUTINE bdy_dyn_frs( kt ) ! Empty routine 180 WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt 199 181 END SUBROUTINE bdy_dyn_frs 200 201 SUBROUTINE bdy_dyn_fla 202 ! No Unstructured open boundaries ==> empty routine 182 SUBROUTINE bdy_dyn_fla ! Empty routine 183 WRITE(*,*) 'bdy_dyn_fla: You should not have seen this print! error?' 203 184 END SUBROUTINE bdy_dyn_fla 204 205 185 #endif 206 186 187 !!====================================================================== 207 188 END MODULE bdydyn -
trunk/NEMO/OPA_SRC/BDY/bdyini.F90
r911 r1125 1 2 MODULE bdyini 3 !!================================================================================= 1 MODULE bdyini 2 !!====================================================================== 4 3 !! *** MODULE bdyini *** 5 !! Initialization of unstructured open boundaries 6 !!================================================================================= 7 #if defined key_bdy || defined key_bdy_tides 8 !!--------------------------------------------------------------------------------- 9 !! 'key_bdy' Unstructured Open Boundary Conditions 10 !!--------------------------------------------------------------------------------- 4 !! Unstructured open boundaries : initialisation 5 !!====================================================================== 6 !! History : 1.0 ! 2005-01 (J. Chanut, A. Sellar) Original code 7 !! - ! 2007-01 (D. Storkey) Update to use IOM module 8 !! - ! 2007-01 (D. Storkey) Tidal forcing 9 !! 3.0 ! 2008-04 (NEMO team) add in the reference version 10 !!---------------------------------------------------------------------- 11 #if defined key_bdy 12 !!---------------------------------------------------------------------- 13 !! 'key_bdy' Unstructured Open Boundary Conditions 14 !!---------------------------------------------------------------------- 11 15 !! bdy_init : Initialization of unstructured open boundaries 12 !!--------------------------------------------------------------------------------- 13 !! * Modules used 16 !!---------------------------------------------------------------------- 14 17 USE oce ! ocean dynamics and tracers variables 15 18 USE dom_oce ! ocean space and time domain … … 24 27 PRIVATE 25 28 26 !! * Routine accessibility 27 PUBLIC bdy_init ! routine called by opa.F90 28 29 !! * Substitutions 30 31 !!--------------------------------------------------------------------------------- 32 !! OPA 9.0 , LODYC-IPSL (2003) 29 PUBLIC bdy_init ! routine called by opa.F90 30 31 !!---------------------------------------------------------------------- 32 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 33 !! $Id: $ 34 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 33 35 !!--------------------------------------------------------------------------------- 34 36 … … 47 49 !! ** Input : bdy_init.nc, input file for unstructured open boundaries 48 50 !! 49 !! History :50 !! OPA 9.0 ! 05-01 (J. Chanut, A. Sellar) Original code51 !! ! 07-01 (D. Storkey) Update to use IOM module.52 !! ! 07-01 (D. Storkey) Tidal forcing.53 51 !!---------------------------------------------------------------------- 54 !! * Local declarations 55 INTEGER :: ji, jj, jk, jgrd, & ! dummy loop indices 56 jb, jr, icount, & 57 icountr, nb_rim, nb_len, nbr_max 52 INTEGER :: ii, ij, ik, igrd, ib, ir ! dummy loop indices 53 INTEGER :: icount, icountr 54 INTEGER :: ib_len, ibr_max 58 55 INTEGER :: iw, ie, is, in 59 56 INTEGER :: inum ! temporary logical unit 60 INTEGER :: & ! temporary integers 61 dummy_id 62 INTEGER, DIMENSION (2) :: kdimsz 63 INTEGER, DIMENSION(jpbdta, jpbgrd) :: & !: Index arrays 64 nbidta, nbjdta, & !: i and j indices of bdy dta 65 nbrdta !: Discrete distance from rim points 66 REAL(wp) :: & 67 efl, wfl, nfl, sfl ! temporary scalars 68 REAL(wp) , DIMENSION(jpidta,jpjdta) :: & 69 tmpmsk ! global domain mask 70 REAL(wp) , DIMENSION(jpbdta,1) :: & 71 ndta ! temporary array 72 CHARACTER(LEN=80),DIMENSION(3) :: bdyfile 73 74 NAMELIST/nambdy/filbdy_mask, filbdy_data_T, filbdy_data_U, filbdy_data_V, & 75 ln_bdy_clim, ln_bdy_vol, ln_bdy_fla, ln_bdy_mask, & 76 nbdy_dta, nb_rimwidth, volbdy 77 57 INTEGER :: id_dummy ! temporary integers 58 INTEGER :: igrd_start, igrd_end ! start and end of loops on igrd 59 INTEGER, DIMENSION (2) :: kdimsz 60 INTEGER, DIMENSION(jpbdta, jpbgrd) :: nbidta, nbjdta ! Index arrays: i and j indices of bdy dta 61 INTEGER, DIMENSION(jpbdta, jpbgrd) :: nbrdta ! Discrete distance from rim points 62 REAL(wp) :: zefl, zwfl, znfl, zsfl ! temporary scalars 63 REAL(wp) , DIMENSION(jpidta,jpjdta) :: zmask ! global domain mask 64 REAL(wp) , DIMENSION(jpbdta,1) :: zdta ! temporary array 65 CHARACTER(LEN=80),DIMENSION(3) :: clfile 66 !! 67 NAMELIST/nambdy/filbdy_mask, filbdy_data_T, filbdy_data_U, filbdy_data_V, & 68 & ln_bdy_tides, ln_bdy_clim, ln_bdy_vol, ln_bdy_mask, & 69 & ln_bdy_dyn_fla, ln_bdy_dyn_frs, ln_bdy_tra_frs, & 70 & nbdy_dta , nb_rimwidth , volbdy 78 71 !!---------------------------------------------------------------------- 79 72 … … 81 74 IF(lwp) WRITE(numout,*) 'bdy_init : initialization of unstructured open boundaries' 82 75 IF(lwp) WRITE(numout,*) '~~~~~~~~' 83 84 IF( jperio /= 0 ) THEN 85 IF(lwp) WRITE(numout,*) 86 IF(lwp) WRITE(numout,*) ' E R R O R : Cyclic or symmetric,', & 87 ' and unstructured open boundary condition are not compatible' 88 IF(lwp) WRITE(numout,*) ' ========== ' 89 IF(lwp) WRITE(numout,*) 90 nstop = nstop + 1 91 END IF 76 ! 77 IF( jperio /= 0 ) CALL ctl_stop( 'Cyclic or symmetric,', & 78 ' and unstructured open boundary condition are not compatible' ) 92 79 93 80 #if defined key_obc 94 IF(lwp) WRITE(numout,*) 95 IF(lwp) WRITE(numout,*) ' E R R O R : Straight open boundaries,', & 96 ' and unstructured open boundaries are not compatible' 97 IF(lwp) WRITE(numout,*) ' ========== ' 98 IF(lwp) WRITE(numout,*) 99 nstop = nstop + 1 81 CALL ctl_stop( 'Straight open boundaries,', & 82 ' and unstructured open boundaries are not compatible' ) 100 83 #endif 101 84 102 85 # if defined key_dynspg_rl 103 IF(lwp) WRITE(numout,*) 104 IF(lwp) WRITE(numout,*) ' E R R O R : Rigid lid,', & 105 ' and unstructured open boundaries are not compatible' 106 IF(lwp) WRITE(numout,*) ' ========== ' 107 IF(lwp) WRITE(numout,*) 108 nstop = nstop + 1 86 CALL ctl_stop( 'Rigid lid,', & 87 ' and unstructured open boundaries are not compatible' ) 109 88 #endif 110 89 111 ! 0.Read namelist parameters90 ! Read namelist parameters 112 91 ! --------------------------- 113 114 92 REWIND( numnam ) 115 93 READ ( numnam, nambdy ) … … 118 96 IF(lwp) WRITE(numout,*) ' nambdy' 119 97 120 IF ((nbdy_dta/=0).AND.(nbdy_dta/=1)) THEN ! Check nbdy_dta value 121 IF(lwp) WRITE(numout,*) 122 IF(lwp) WRITE(numout,*) ' E R R O R : nbdy_dta =',nbdy_dta, & 123 'but it should have been 0 or 1' 124 IF(lwp) WRITE(numout,*) ' ========== ' 125 IF(lwp) WRITE(numout,*) 126 nstop = nstop + 1 127 ELSE 128 IF(lwp) WRITE(numout,*) ' ' 129 IF(lwp) WRITE(numout,*) ' data in file (=1) or nbdy_dta = ', nbdy_dta 130 IF(lwp) WRITE(numout,*) ' initial state used (=0)' 131 END IF 98 ! Check nbdy_dta value 99 IF(lwp) WRITE(numout,*) 'nbdy_dta =', nbdy_dta 100 IF(lwp) WRITE(numout,*) ' ' 101 SELECT CASE( nbdy_dta ) 102 CASE( 0 ) 103 IF(lwp) WRITE(numout,*) ' initial state used for bdy data' 104 CASE( 1 ) 105 IF(lwp) WRITE(numout,*) ' boundary data taken from file' 106 CASE DEFAULT 107 CALL ctl_stop( 'nbdy_dta must be 0 or 1' ) 108 END SELECT 132 109 133 110 IF(lwp) WRITE(numout,*) ' ' 134 111 IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nb_rimwidth = ', nb_rimwidth 135 112 113 IF(lwp) WRITE(numout,*) ' ' 114 IF(lwp) WRITE(numout,*) ' volbdy = ', volbdy 115 136 116 IF (ln_bdy_vol) THEN 137 IF (volbdy==1) THEN ! Check volbdy value 138 IF(lwp) WRITE(numout,*) ' ' 139 IF(lwp) WRITE(numout,*) ' volbdy = ', volbdy 117 SELECT CASE ( volbdy ) ! Check volbdy value 118 CASE( 1 ) 140 119 IF(lwp) WRITE(numout,*) ' The total volume will be constant' 141 142 ELSEIF (volbdy==0) THEN 143 IF(lwp) WRITE(numout,*) ' ' 144 IF(lwp) WRITE(numout,*) ' volbdy = ', volbdy 145 IF(lwp) WRITE(numout,*) ' The total volume will vary according to & 146 &the surface E-P flux' 147 ELSE 148 IF(lwp) WRITE(numout,*) ' ' 149 IF(lwp) WRITE(numout,*) ' E R R O R : volbdy =',volbdy, & 150 'but it should have been 0 or 1' 151 IF(lwp) WRITE(numout,*) ' ========== ' 152 IF(lwp) WRITE(numout,*) 153 nstop = nstop + 1 154 END IF 120 CASE( 0 ) 121 IF(lwp) WRITE(numout,*) ' The total volume will vary according' 122 IF(lwp) WRITE(numout,*) ' to the surface E-P flux' 123 CASE DEFAULT 124 CALL ctl_stop( 'volbdy must be 0 or 1' ) 125 END SELECT 155 126 ELSE 156 IF(lwp) WRITE(numout,*) ' '157 127 IF(lwp) WRITE(numout,*) 'No volume correction with unstructured open boundaries' 158 128 IF(lwp) WRITE(numout,*) ' ' 159 129 ENDIF 160 130 161 IF (ln_bdy_fla) THEN 162 IF(lwp) WRITE(numout,*) ' ' 163 IF(lwp) WRITE(numout,*) 'Flather bc with unstructured open boundaries' 164 IF(lwp) WRITE(numout,*) ' ' 165 ELSE 166 IF(lwp) WRITE(numout,*) ' ' 167 IF(lwp) WRITE(numout,*) 'NO Flather bc with unstructured open boundaries' 168 IF(lwp) WRITE(numout,*) ' ' 169 ENDIF 170 171 ! 0.5 Read tides namelist 131 IF (ln_bdy_tides) THEN 132 IF(lwp) WRITE(numout,*) ' ' 133 IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing at unstructured open boundaries' 134 IF(lwp) WRITE(numout,*) ' ' 135 ENDIF 136 137 IF (ln_bdy_dyn_fla) THEN 138 IF(lwp) WRITE(numout,*) ' ' 139 IF(lwp) WRITE(numout,*) 'Flather condition on U, V at unstructured open boundaries' 140 IF(lwp) WRITE(numout,*) ' ' 141 ENDIF 142 143 IF (ln_bdy_dyn_frs) THEN 144 IF(lwp) WRITE(numout,*) ' ' 145 IF(lwp) WRITE(numout,*) 'FRS condition on U and V at unstructured open boundaries' 146 IF(lwp) WRITE(numout,*) ' ' 147 ENDIF 148 149 IF (ln_bdy_tra_frs) THEN 150 IF(lwp) WRITE(numout,*) ' ' 151 IF(lwp) WRITE(numout,*) 'FRS condition on T & S fields at unstructured open boundaries' 152 IF(lwp) WRITE(numout,*) ' ' 153 ENDIF 154 155 ! Read tides namelist 172 156 ! ------------------------ 173 174 IF ( lk_bdy_tides ) CALL tide_init 175 176 ! 1. Read arrays defining unstructured open boundaries 177 ! ---------------------------------------------------- 178 179 ! 1.1 Read global 2D mask at T-points: bdytmask 180 ! ********************************************* 181 ! bdytmask=1 on the computational domain AND on open boundaries 182 ! =0 elsewhere 157 IF( ln_bdy_tides ) CALL tide_init 158 159 ! Read arrays defining unstructured open boundaries 160 ! ------------------------------------------------- 161 162 ! Read global 2D mask at T-points: bdytmask 163 ! ***************************************** 164 ! bdytmask = 1 on the computational domain AND on open boundaries 165 ! = 0 elsewhere 183 166 184 167 IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN 185 tmpmsk(: , :) = 0.e0186 tmpmsk(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0168 zmask( : ,:) = 0.e0 169 zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0 187 170 ELSE IF ( ln_bdy_mask ) THEN 188 CALL iom_open( filbdy_mask, inum )189 CALL iom_get ( inum, jpdom_data, 'bdy_msk', tmpmsk(:,:) )190 CALL iom_close( inum )171 CALL iom_open( filbdy_mask, inum ) 172 CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) ) 173 CALL iom_close( inum ) 191 174 ELSE 192 tmpmsk(:,:) = 1.0175 zmask(:,:) = 1.e0 193 176 ENDIF 194 177 195 178 ! Save mask over local domain 196 DO jj = 1, nlcj197 DO ji = 1, nlci198 bdytmask(ji,jj) = tmpmsk( mig(ji), mjg(jj))199 END DO179 DO ij = 1, nlcj 180 DO ii = 1, nlci 181 bdytmask(ii,ij) = zmask( mig(ii), mjg(ij) ) 182 END DO 200 183 END DO 201 184 202 185 ! Derive mask on U and V grid from mask on T grid 203 bdyumask(:,:) =0.e0204 bdyvmask(:,:) =0.e0205 DO jj=1, jpjm1206 DO ji=1, jpim1207 bdyumask(ji,jj)=bdytmask(ji,jj)*bdytmask(ji+1, jj )208 bdyvmask(ji,jj)=bdytmask(ji,jj)*bdytmask(ji ,jj+1)209 END DO186 bdyumask(:,:) = 0.e0 187 bdyvmask(:,:) = 0.e0 188 DO ij=1, jpjm1 189 DO ii=1, jpim1 190 bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij ) 191 bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii ,ij+1) 192 END DO 210 193 END DO 211 194 … … 214 197 CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) 215 198 216 ! 1.2Read discrete distance and mapping indices217 ! ****************************************** ****218 nbidta(:,:)=0.219 nbjdta(:,:)=0.220 nbrdta(:,:)=0.199 ! Read discrete distance and mapping indices 200 ! ****************************************** 201 nbidta(:,:) = 0.e0 202 nbjdta(:,:) = 0.e0 203 nbrdta(:,:) = 0.e0 221 204 222 205 IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN 223 224 icount = 0 225 ! Define west boundary (from ji=2 to ji=1+nb_rimwidth): 226 DO jr=1,nb_rimwidth 227 DO jj=3,jpjglo-2 228 icount=icount+1 229 nbidta(icount,:) = jr + 1 + (jpizoom-1) 230 nbjdta(icount,:) = jj + (jpjzoom-1) 231 nbrdta(icount,:) = jr 232 END DO 233 END DO 234 235 ! Define east boundary (from ji=jpiglo-1 to ji=jpiglo-nb_rimwidth): 236 DO jr=1,nb_rimwidth 237 DO jj=3,jpjglo-2 238 icount=icount+1 239 nbidta(icount,:) = jpiglo-jr + (jpizoom-1) 240 nbidta(icount,2) = jpiglo-jr-1 + (jpizoom-1) ! special case for u points 241 nbjdta(icount,:) = jj + (jpjzoom-1) 242 nbrdta(icount,:) = jr 243 END DO 244 END DO 206 icount = 0 207 ! Define west boundary (from ii=2 to ii=1+nb_rimwidth): 208 DO ir = 1, nb_rimwidth 209 DO ij = 3, jpjglo-2 210 icount=icount+1 211 nbidta(icount,:) = ir + 1 + (jpizoom-1) 212 nbjdta(icount,:) = ij + (jpjzoom-1) 213 nbrdta(icount,:) = ir 214 END DO 215 END DO 216 217 ! Define east boundary (from ii=jpiglo-1 to ii=jpiglo-nb_rimwidth): 218 DO ir=1,nb_rimwidth 219 DO ij=3,jpjglo-2 220 icount=icount+1 221 nbidta(icount,:) = jpiglo-ir + (jpizoom-1) 222 nbidta(icount,2) = jpiglo-ir-1 + (jpizoom-1) ! special case for u points 223 nbjdta(icount,:) = ij + (jpjzoom-1) 224 nbrdta(icount,:) = ir 225 END DO 226 END DO 245 227 246 ELSE 247 248 ! Read indices and distances in unstructured boundary data files 249 250 IF ( lk_bdy ) THEN 251 bdyfile(1) = filbdy_data_T 252 bdyfile(2) = filbdy_data_U 253 bdyfile(3) = filbdy_data_V 254 ELSE 255 ! In this case we have tides only at the boundaries 256 ! so read index arrays from tides files for first tidal component 257 bdyfile(1) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_T.nc' 258 bdyfile(2) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_U.nc' 259 bdyfile(3) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_V.nc' 260 ENDIF 261 262 DO jgrd = 1,3 263 CALL iom_open( bdyfile(jgrd), inum ) 264 dummy_id = iom_varid( inum, 'nbidta', kdimsz=kdimsz ) 265 WRITE(numout,*) 'kdimsz : ',kdimsz 266 nb_len = kdimsz(1) 267 IF (nb_len > jpbdta) THEN 268 IF(lwp) WRITE(numout,*) 269 IF(lwp) WRITE(numout,*) ' E R R O R : jpbdta is too small:' 270 IF(lwp) WRITE(numout,*) ' ========== Boundary array length in file is ', nb_len 271 IF(lwp) WRITE(numout,*) ' But jpbdta is ', jpbdta 272 IF(lwp) WRITE(numout,*) ' File : ', bdyfile(jgrd) 273 IF(lwp) WRITE(numout,*) 274 nstop = nstop + 1 275 ENDIF 276 CALL iom_get ( inum, jpdom_unknown, 'nbidta', ndta(1:nb_len,:) ) 277 DO ji=1,nb_len 278 nbidta(ji,jgrd) = INT( ndta(ji,1) ) 279 ENDDO 280 CALL iom_get ( inum, jpdom_unknown, 'nbjdta', ndta(1:nb_len,:) ) 281 DO ji=1,nb_len 282 nbjdta(ji,jgrd) = INT( ndta(ji,1) ) 283 ENDDO 284 CALL iom_get ( inum, jpdom_unknown, 'nbrdta', ndta(1:nb_len,:) ) 285 DO ji=1,nb_len 286 nbrdta(ji,jgrd) = INT( ndta(ji,1) ) 287 ENDDO 288 CALL iom_close( inum ) 289 290 ! Check that rimwidth in file is big enough: 291 nbr_max = MAXVAL(nbrdta(:,jgrd)) 292 IF (nbr_max < nb_rimwidth) THEN 293 IF(lwp) WRITE(numout,*) 294 IF(lwp) WRITE(numout,*) ' E R R O R : Maximum rimwidth in file is ', nbr_max 295 IF(lwp) WRITE(numout,*) ' ========== but nb_rimwidth is ', nb_rimwidth 296 IF(lwp) WRITE(numout,*) 297 nstop = nstop + 1 298 ELSE 299 IF(lwp) WRITE(numout,*) 300 IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', nbr_max 301 IF(lwp) WRITE(numout,*) 302 END IF 303 304 ENDDO 305 306 END IF 307 308 ! 1.3 Dispatch mapping indices and discrete distances on each processor 309 ! ********************************************************************* 228 ELSE ! Read indices and distances in unstructured boundary data files 229 230 IF( ln_bdy_tides ) THEN 231 ! Read tides input files for preference in case there are 232 ! no bdydata files. 233 clfile(1) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_T.nc' 234 clfile(2) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_U.nc' 235 clfile(3) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_V.nc' 236 ELSE 237 clfile(1) = filbdy_data_T 238 clfile(2) = filbdy_data_U 239 clfile(3) = filbdy_data_V 240 ENDIF 241 242 ! how many files are we to read in? 243 igrd_start = 1 244 igrd_end = 3 245 IF(.NOT. ln_bdy_tides ) THEN 246 IF(.NOT. (ln_bdy_dyn_fla) .AND..NOT. (ln_bdy_tra_frs)) THEN 247 ! No T-grid file. 248 igrd_start = 2 249 ELSEIF ( .NOT. ln_bdy_dyn_frs .AND..NOT. ln_bdy_dyn_fla ) THEN 250 ! No U-grid or V-grid file. 251 igrd_end = 1 252 ENDIF 253 ENDIF 254 255 DO igrd = igrd_start, igrd_end 256 CALL iom_open( clfile(igrd), inum ) 257 id_dummy = iom_varid( inum, 'nbidta', kdimsz=kdimsz ) 258 WRITE(numout,*) 'kdimsz : ',kdimsz 259 ib_len = kdimsz(1) 260 IF( ib_len > jpbdta) CALL ctl_stop( & 261 'Boundary data array in file too long.', & 262 'File :', TRIM(clfile(igrd)), & 263 'increase parameter jpbdta.' ) 264 265 CALL iom_get( inum, jpdom_unknown, 'nbidta', zdta(1:ib_len,:) ) 266 DO ii = 1,ib_len 267 nbidta(ii,igrd) = INT( zdta(ii,1) ) 268 END DO 269 CALL iom_get( inum, jpdom_unknown, 'nbjdta', zdta(1:ib_len,:) ) 270 DO ii = 1,ib_len 271 nbjdta(ii,igrd) = INT( zdta(ii,1) ) 272 END DO 273 CALL iom_get ( inum, jpdom_unknown, 'nbrdta', zdta(1:ib_len,:) ) 274 DO ii = 1,ib_len 275 nbrdta(ii,igrd) = INT( zdta(ii,1) ) 276 END DO 277 CALL iom_close( inum ) 278 279 ! Check that rimwidth in file is big enough: 280 ibr_max = MAXVAL( nbrdta(:,igrd) ) 281 IF(lwp) WRITE(numout,*) 282 IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max 283 IF(lwp) WRITE(numout,*) ' nb_rimwidth from namelist is ', nb_rimwidth 284 IF (ibr_max < nb_rimwidth) CALL ctl_stop( & 285 'nb_rimwidth is larger than maximum rimwidth in file' ) 286 ! 287 END DO 288 ! 289 ENDIF 290 291 ! Dispatch mapping indices and discrete distances on each processor 292 ! ***************************************************************** 310 293 311 iw = mig(1) +1 ! if monotasking and no zoom, iw=2312 ie = mig(1) + nlci-1 -1 ! if monotasking and no zoom, ie=jpim1313 is = mjg(1) +1 ! if monotasking and no zoom, is=2314 in = mjg(1) + nlcj-1 -1 ! if monotasking and no zoom, in=jpjm1315 316 DO jgrd = 1, jpbgrd294 iw = mig(1) + 1 ! if monotasking and no zoom, iw=2 295 ie = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1 296 is = mjg(1) + 1 ! if monotasking and no zoom, is=2 297 in = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1 298 299 DO igrd = igrd_start, igrd_end 317 300 icount = 0 318 301 icountr = 0 319 DO nb_rim=1, nb_rimwidth320 DO jb = 1, jpbdta321 ! check if point is in local domain and equals nb_rim322 IF ( (nbidta(jb,jgrd) >= iw ).AND. &323 (nbidta(jb,jgrd) <= ie ).AND. &324 (nbjdta(jb,jgrd) >= is ).AND. &325 (nbjdta(jb,jgrd) <= in ).AND.&326 (nbrdta(jb,jgrd) == nb_rim ) ) THEN327 328 icount = icount + 1329 330 IF (nb_rim==1) icountr = icountr+1331 332 IF (icount > jpbdim) THEN333 IF(lwp) WRITE(numout,*) 'bdy_ini: jpbdim too small'334 nstop = nstop + 1335 ELSE336 nbi(icount, jgrd) = nbidta(jb,jgrd)- mig(1)+1337 nbj(icount, jgrd) = nbjdta(jb,jgrd)- mjg(1)+1338 nbr(icount, jgrd) = nbrdta(jb,jgrd)339 nbmap(icount,jgrd) = jb340 ENDIF341 ENDIF342 END DO343 END DO344 nblenrim(jgrd) = icountr !: length of rim boundary data on each proc345 nblen (jgrd) = icount !: length of boundary data on each proc302 nblen(igrd) = 0 303 nblenrim(igrd) = 0 304 nblendta(igrd) = 0 305 DO ir=1, nb_rimwidth 306 DO ib = 1, jpbdta 307 ! check if point is in local domain and equals ir 308 IF( nbidta(ib,igrd) >= iw .AND. nbidta(ib,igrd) <= ie .AND. & 309 & nbjdta(ib,igrd) >= is .AND. nbjdta(ib,igrd) <= in .AND. & 310 & nbrdta(ib,igrd) == ir ) THEN 311 ! 312 icount = icount + 1 313 ! 314 IF( ir == 1 ) icountr = icountr+1 315 IF (icount > jpbdim) THEN 316 IF(lwp) WRITE(numout,*) 'bdy_ini: jpbdim too small' 317 nstop = nstop + 1 318 ELSE 319 nbi(icount, igrd) = nbidta(ib,igrd)- mig(1)+1 320 nbj(icount, igrd) = nbjdta(ib,igrd)- mjg(1)+1 321 nbr(icount, igrd) = nbrdta(ib,igrd) 322 nbmap(icount,igrd) = ib 323 ENDIF 324 ENDIF 325 END DO 326 END DO 327 nblenrim(igrd) = icountr !: length of rim boundary data on each proc 328 nblen (igrd) = icount !: length of boundary data on each proc 346 329 END DO 347 330 348 ! 2. Compute rim weights 349 ! ---------------------- 350 351 DO jgrd = 1, jpbgrd 352 DO jb = 1, nblen(jgrd) 353 ! tanh formulation 354 nbw(jb,jgrd) = 1.-TANH((FLOAT(nbr(jb,jgrd)-1))/2.) 355 ! quadratic 356 ! nbw(jb,jgrd) = (FLOAT(nb_rimwidth+1-nbr(jb,jgrd))/FLOAT(nb_rimwidth))**2 357 ! linear 358 ! nbw(jb,jgrd) = FLOAT(nb_rimwidth+1-nbr(jb,jgrd))/FLOAT(nb_rimwidth) 359 END DO 331 ! Compute rim weights 332 ! ------------------- 333 DO igrd = igrd_start, igrd_end 334 DO ib = 1, nblen(igrd) 335 ! tanh formulation 336 nbw(ib,igrd) = 1.- TANH( FLOAT( nbr(ib,igrd) - 1 ) *0.5 ) 337 ! quadratic 338 ! nbw(ib,igrd) = (FLOAT(nb_rimwidth+1-nbr(ib,igrd))/FLOAT(nb_rimwidth))**2 339 ! linear 340 ! nbw(ib,igrd) = FLOAT(nb_rimwidth+1-nbr(ib,igrd))/FLOAT(nb_rimwidth) 341 END DO 360 342 END DO 361 343 362 ! 3. Mask corrections 363 ! ------------------- 364 365 DO jk=1, jpkm1 366 DO jj=1, jpj 367 DO ji=1, jpi 368 tmask(ji,jj,jk)=tmask(ji,jj,jk)*bdytmask(ji,jj) 369 umask(ji,jj,jk)=umask(ji,jj,jk)*bdyumask(ji,jj) 370 vmask(ji,jj,jk)=vmask(ji,jj,jk)*bdyvmask(ji,jj) 371 bmask(ji,jj)=bmask(ji,jj)*bdytmask(ji,jj) 372 END DO 373 END DO 374 END DO 375 376 ! I am not sure that it is useful: 377 DO jk=1, jpkm1 378 DO jj=2, jpjm1 379 DO ji=2, jpim1 380 fmask(ji,jj,jk) = fmask(ji,jj,jk) & 381 & * bdytmask(ji, jj ) * bdytmask(ji+1, jj ) & 382 & * bdytmask(ji,jj+1) * bdytmask(ji+1,jj+1) 383 END DO 384 END DO 385 END DO 386 387 tmask_i(:,:) = tmask(:,:,1)*tmask_i(:,:) 388 389 bdytmask(:,:)=tmask(:,:,1) 344 ! Mask corrections 345 ! ---------------- 346 DO ik = 1, jpkm1 347 DO ij = 1, jpj 348 DO ii = 1, jpi 349 tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij) 350 umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij) 351 vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij) 352 bmask(ii,ij) = bmask(ii,ij) * bdytmask(ii,ij) 353 END DO 354 END DO 355 END DO 356 357 DO ik = 1, jpkm1 358 DO ij = 2, jpjm1 359 DO ii = 2, jpim1 360 fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij ) * bdytmask(ii+1,ij ) & 361 & * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1) 362 END DO 363 END DO 364 END DO 365 366 tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:) 367 bdytmask(:,:) = tmask(:,:,1) 390 368 391 369 ! bdy masks and bmask are now set to zero on boundary points: 392 393 jgrd=1 ! In the free surface case, bmask is at T-points 394 DO jb=1, nblenrim(jgrd) 395 bmask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0 396 END DO 397 398 jgrd=1 399 DO jb=1, nblenrim(jgrd) 400 bdytmask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0 401 END DO 402 403 jgrd=2 404 DO jb=1, nblenrim(jgrd) 405 bdyumask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0 406 END DO 407 408 jgrd=3 409 DO jb=1, nblenrim(jgrd) 410 bdyvmask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0 411 END DO 412 413 ! Lateral boundary conditions 414 415 CALL lbc_lnk( fmask, 'F', 1. ) 370 igrd = 1 ! In the free surface case, bmask is at T-points 371 DO ib = 1, nblenrim(igrd) 372 bmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 373 END DO 374 ! 375 igrd = 1 376 DO ib = 1, nblenrim(igrd) 377 bdytmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 378 END DO 379 ! 380 igrd = 2 381 DO ib = 1, nblenrim(igrd) 382 bdyumask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 383 END DO 384 ! 385 igrd = 3 386 DO ib = 1, nblenrim(igrd) 387 bdyvmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 388 END DO 389 390 ! Lateral boundary conditions 391 CALL lbc_lnk( fmask , 'F', 1. ) 416 392 CALL lbc_lnk( bdytmask(:,:), 'T', 1. ) 417 393 CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) 418 394 CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) 419 395 420 IF ((ln_bdy_vol).OR.(ln_bdy_fla)) THEN 421 422 ! 4 Indices and directions of rim velocity components 423 ! --------------------------------------------------- 424 !flagu = -1 : u component is normal to the dynamical boundary 425 ! but its direction is outward 426 ! 427 !flagu = 0 : u is tangential 428 ! 429 !flagu = 1 : u is normal to the boundary 430 ! and is direction is inward 396 IF( ln_bdy_vol .OR. ln_bdy_dyn_fla ) THEN ! Indices and directions of rim velocity components 397 ! 398 !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward 399 !flagu = 0 : u is tangential 400 !flagu = 1 : u is normal to the boundary and is direction is inward 401 icount = 0 402 flagu(:) = 0.e0 431 403 432 icount = 0 433 434 flagu(:)=0.e0 404 igrd = 2 ! u-component 405 DO ib = 1, nblenrim(igrd) 406 zefl=bdytmask(nbi(ib,igrd) , nbj(ib,igrd)) 407 zwfl=bdytmask(nbi(ib,igrd)+1, nbj(ib,igrd)) 408 IF( zefl + zwfl ==2 ) THEN 409 icount = icount +1 410 ELSE 411 flagu(ib)=-zefl+zwfl 412 ENDIF 413 END DO 414 415 !flagv = -1 : u component is normal to the dynamical boundary but its direction is outward 416 !flagv = 0 : u is tangential 417 !flagv = 1 : u is normal to the boundary and is direction is inward 418 flagv(:) = 0.e0 419 420 igrd = 3 ! v-component 421 DO ib = 1, nblenrim(igrd) 422 znfl = bdytmask(nbi(ib,igrd), nbj(ib,igrd)) 423 zsfl = bdytmask(nbi(ib,igrd), nbj(ib,igrd)+1) 424 IF( znfl + zsfl ==2 ) THEN 425 icount = icount + 1 426 ELSE 427 flagv(ib) = -znfl + zsfl 428 END IF 429 END DO 435 430 436 jgrd=2 ! u-component 437 DO jb=1, nblenrim(jgrd) 438 efl=bdytmask(nbi(jb,jgrd) , nbj(jb,jgrd)) 439 wfl=bdytmask(nbi(jb,jgrd)+1, nbj(jb,jgrd)) 440 IF ((efl+wfl)==2) THEN 441 icount = icount +1 442 ELSE 443 flagu(jb)=-efl+wfl 444 END IF 445 END DO 446 447 !flagv = -1 : u component is normal to the dynamical boundary 448 ! but its direction is outward 449 ! 450 !flagv = 0 : u is tangential 451 ! 452 !flagv = 1 : u is normal to the boundary 453 ! and is direction is inward 454 455 flagv(:)=0.e0 456 457 jgrd=3 ! v-component 458 DO jb=1, nblenrim(jgrd) 459 nfl = bdytmask(nbi(jb,jgrd), nbj(jb,jgrd)) 460 sfl = bdytmask(nbi(jb,jgrd), nbj(jb,jgrd)+1) 461 IF ((nfl+sfl)==2) THEN 462 icount = icount +1 463 ELSE 464 flagv(jb)=-nfl+sfl 465 END IF 466 END DO 467 468 IF( icount /= 0 ) THEN 469 IF(lwp) WRITE(numout,*) 470 IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,', & 471 ' are not boundary points. Check nbi, nbj, indices.' 472 IF(lwp) WRITE(numout,*) ' ========== ' 473 IF(lwp) WRITE(numout,*) 474 nstop = nstop + 1 475 END IF 431 IF( icount /= 0 ) THEN 432 IF(lwp) WRITE(numout,*) 433 IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,', & 434 ' are not boundary points. Check nbi, nbj, indices.' 435 IF(lwp) WRITE(numout,*) ' ========== ' 436 IF(lwp) WRITE(numout,*) 437 nstop = nstop + 1 438 ENDIF 476 439 477 END 478 479 ! 5Compute total lateral surface for volume correction:480 ! ---------------------------------------------------- --440 ENDIF 441 442 ! Compute total lateral surface for volume correction: 443 ! ---------------------------------------------------- 481 444 482 445 bdysurftot = 0.e0 483 484 IF (ln_bdy_vol) THEN 485 jgrd=2 ! Lateral surface at U-points 486 DO jb=1, nblenrim(jgrd) 487 bdysurftot = bdysurftot + & 488 hu(nbi(jb,jgrd), nbj(jb,jgrd)) * e2u(nbi(jb,jgrd), nbj(jb,jgrd)) & 489 * ABS(flagu(jb))*tmask_i(nbi(jb,jgrd) , nbj(jb,jgrd)) & 490 *tmask_i(nbi(jb,jgrd)+1, nbj(jb,jgrd)) 491 END DO 492 493 jgrd=3 ! Add lateral surface at V-points 494 DO jb=1, nblenrim(jgrd) 495 bdysurftot = bdysurftot + & 496 hv(nbi(jb,jgrd), nbj(jb,jgrd)) * e1v(nbi(jb,jgrd), nbj(jb,jgrd)) & 497 * ABS(flagv(jb))*tmask_i(nbi(jb,jgrd), nbj(jb,jgrd)) & 498 *tmask_i(nbi(jb,jgrd), nbj(jb,jgrd)+1) 499 END DO 500 501 IF( lk_mpp ) CALL mpp_sum( bdysurftot ) ! sum over the global domain 446 IF( ln_bdy_vol ) THEN 447 igrd = 2 ! Lateral surface at U-points 448 DO ib = 1, nblenrim(igrd) 449 bdysurftot = bdysurftot + hu (nbi(ib,igrd) ,nbj(ib,igrd)) & 450 & * e2u (nbi(ib,igrd) ,nbj(ib,igrd)) * ABS( flagu(ib) ) & 451 & * tmask_i(nbi(ib,igrd) ,nbj(ib,igrd)) & 452 & * tmask_i(nbi(ib,igrd)+1,nbj(ib,igrd)) 453 END DO 454 455 igrd=3 ! Add lateral surface at V-points 456 DO ib = 1, nblenrim(igrd) 457 bdysurftot = bdysurftot + hv (nbi(ib,igrd),nbj(ib,igrd) ) & 458 & * e1v (nbi(ib,igrd),nbj(ib,igrd) ) * ABS( flagv(ib) ) & 459 & * tmask_i(nbi(ib,igrd),nbj(ib,igrd) ) & 460 & * tmask_i(nbi(ib,igrd),nbj(ib,igrd)+1) 461 END DO 462 463 IF( lk_mpp ) CALL mpp_sum( bdysurftot ) ! sum over the global domain 502 464 END IF 503 465 504 505 ! 6. Initialise bdy data arrays 506 ! ----------------------------- 507 466 ! Initialise bdy data arrays 467 ! -------------------------- 508 468 tbdy(:,:) = 0.e0 509 469 sbdy(:,:) = 0.e0 … … 514 474 vbtbdy(:) = 0.e0 515 475 516 ! 7. Read in tidal constituents and adjust for model start time 517 ! ------------------------------------------------------------- 518 519 IF ( lk_bdy_tides ) CALL tide_data 520 476 ! Read in tidal constituents and adjust for model start time 477 ! ---------------------------------------------------------- 478 IF( ln_bdy_tides ) CALL tide_data 479 ! 521 480 END SUBROUTINE bdy_init 522 481 -
trunk/NEMO/OPA_SRC/BDY/bdytides.F90
r911 r1125 1 1 MODULE bdytides 2 !!====================================================================== ===========2 !!====================================================================== 3 3 !! *** MODULE bdytides *** 4 4 !! Ocean dynamics: Tidal forcing at open boundaries 5 !!================================================================================= 6 #if defined key_bdy_tides 7 !!--------------------------------------------------------------------------------- 5 !!====================================================================== 6 !! History : 2.0 ! 2007-01 (D.Storkey) Original code 7 !! 2.3 ! 2008-01 (J.Holt) Add date correction. Origins POLCOMS v6.3 2007 8 !! 3.0 ! 2008-04 (NEMO team) add in the reference version 9 !!---------------------------------------------------------------------- 10 #if defined key_bdy 11 !!---------------------------------------------------------------------- 12 !! 'key_bdy' Unstructured Open Boundary Condition 13 !!---------------------------------------------------------------------- 8 14 !! PUBLIC 9 !! tide_init: read of namelist10 !! tide_data: read in and initialisation of tidal constituents at boundary11 !! tide_update: calculation of tidal forcing at each timestep15 !! tide_init : read of namelist 16 !! tide_data : read in and initialisation of tidal constituents at boundary 17 !! tide_update : calculation of tidal forcing at each timestep 12 18 !! PRIVATE 13 !! uvset :\ 14 !! vday : | Routines to correct tidal harmonics forcing for 15 !! shpen : | start time of integration 16 !! ufset : | 17 !! vset :/ 18 !!--------------------------------------------------------------------------------- 19 20 !!--------------------------------------------------------------------------------- 21 !! * Modules used 19 !! uvset :\ 20 !! vday : | Routines to correct tidal harmonics forcing for 21 !! shpen : | start time of integration 22 !! ufset : | 23 !! vset :/ 24 !!---------------------------------------------------------------------- 22 25 USE oce ! ocean dynamics and tracers 23 26 USE dom_oce ! ocean space and time domain … … 29 32 USE bdy_oce ! ocean open boundary conditions 30 33 USE daymod ! calendar 34 31 35 IMPLICIT NONE 32 36 PRIVATE 33 37 34 !! * Accessibility 35 PUBLIC tide_init ! routine called in bdyini 36 PUBLIC tide_data ! routine called in bdyini 37 PUBLIC tide_update ! routine called in bdydyn 38 39 !! * Module variables 40 LOGICAL, PUBLIC, PARAMETER :: lk_bdy_tides = .TRUE. !: tidal forcing at boundaries. 41 42 CHARACTER(len=80), PUBLIC :: & 43 filtide !: Filename root for tidal input files 44 45 INTEGER, PUBLIC, PARAMETER :: ntide_max = 15 ! Max number of tidal contituents 46 INTEGER, PUBLIC :: ntide ! Actual number of tidal constituents 47 48 CHARACTER(len=4), PUBLIC, DIMENSION(ntide_max) :: & 49 tide_cpt !: Names of tidal components used. 50 51 logical, PUBLIC :: ln_tide_date ! if true correct tide phases and amplitude for model start date 52 53 REAL(wp), PUBLIC, DIMENSION(ntide_max) :: tide_speed ! Phase speed of tidal constituent (deg/hr) 54 INTEGER, DIMENSION(ntide_max) :: indx 55 REAL(wp), DIMENSION(jpbdim,ntide_max) :: & 56 ssh1, ssh2, & !: Tidal constituents : SSH 57 u1 , u2 , & !: Tidal constituents : U 58 v1 , v2 !: Tidal constituents : V 38 PUBLIC tide_init ! routine called in bdyini 39 PUBLIC tide_data ! routine called in bdyini 40 PUBLIC tide_update ! routine called in bdydyn 41 42 LOGICAL, PUBLIC :: ln_tide_date !: =T correct tide phases and amplitude for model start date 43 44 INTEGER, PARAMETER :: jptides_max = 15 !: Max number of tidal contituents 45 INTEGER :: ntide !: Actual number of tidal constituents 46 47 CHARACTER(len=80), PUBLIC :: filtide !: Filename root for tidal input files 48 CHARACTER(len= 4), PUBLIC, DIMENSION(jptides_max) :: tide_cpt !: Names of tidal components used. 49 50 INTEGER , DIMENSION(jptides_max) :: nindx !: ??? 51 REAL(wp), DIMENSION(jptides_max) :: tide_speed !: Phase speed of tidal constituent (deg/hr) 59 52 53 REAL(wp), DIMENSION(jpbdim,jptides_max) :: ssh1, ssh2 !: Tidal constituents : SSH 54 REAL(wp), DIMENSION(jpbdim,jptides_max) :: u1 , u2 !: Tidal constituents : U 55 REAL(wp), DIMENSION(jpbdim,jptides_max) :: v1 , v2 !: Tidal constituents : V 60 56 61 !!--------------------------------------------------------------------------------- 57 !!---------------------------------------------------------------------- 58 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 59 !! $Id: $ 60 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 61 !!---------------------------------------------------------------------- 62 62 63 63 CONTAINS 64 64 65 65 SUBROUTINE tide_init 66 !!---------------------------------------------------------------------- --------67 !! SUBROUTINE tide_init68 !! ***********************66 !!---------------------------------------------------------------------- 67 !! *** SUBROUTINE tide_init *** 68 !! 69 69 !! ** Purpose : - Read in namelist for tides 70 70 !! 71 !! History : 72 !! NEMO v2.0 ! 07-01 (D.Storkey) Original 73 !!------------------------------------------------------------------------------ 74 !! * Local declarations 75 INTEGER :: jtide ! dummy loop index 76 ! different from nblendta!) 77 78 !!------------------------------------------------------------------------------ 79 !! OPA 9.0, LODYC-IPSL (2007) 80 !!------------------------------------------------------------------------------ 81 82 NAMELIST/namtide/filtide, tide_cpt, tide_speed, ln_tide_date 83 71 !!---------------------------------------------------------------------- 72 INTEGER :: itide ! dummy loop index 73 !! 74 NAMELIST/namtide/ln_tide_date, filtide, tide_cpt, tide_speed 84 75 !!---------------------------------------------------------------------- 85 76 … … 88 79 IF(lwp) WRITE(numout,*) '~~~~~~~~~' 89 80 90 ! 0. Read namelist parameters 91 ! --------------------------- 92 93 do jtide = 1, ntide_max 94 tide_cpt(jtide) = '' 95 enddo 96 97 REWIND( numnam ) 81 ! Namelist namtide : tidal harmonic forcing at open boundaries 82 ln_tide_date = .false. 83 filtide(:) = '' 84 tide_speed(:) = 0.0 85 tide_cpt(:) = '' 86 REWIND( numnam ) ! Read namelist parameters 98 87 READ ( numnam, namtide ) 99 100 ! control prints 101 IF(lwp) WRITE(numout,*) ' namtide' 102 103 ! Count number of components specified: 104 ntide=ntide_max 105 do jtide = 1, ntide_max 106 if ( tide_cpt(jtide) == '' ) then 107 ntide = jtide-1 108 exit 109 endif 110 111 enddo 112 ! find constients in standard list 113 do jtide = 1, ntide 114 indx(jtide) = 0 115 if ( TRIM(tide_cpt(jtide)) == 'Q1' ) indx(jtide) = 1 116 if ( TRIM(tide_cpt(jtide)) == 'O1' ) indx(jtide) = 2 117 if ( TRIM(tide_cpt(jtide)) == 'P1' ) indx(jtide) = 3 118 if ( TRIM(tide_cpt(jtide)) == 'S1' ) indx(jtide) = 4 119 if ( TRIM(tide_cpt(jtide)) == 'K1' ) indx(jtide) = 5 120 if ( TRIM(tide_cpt(jtide)) == '2N2' ) indx(jtide) = 6 121 if ( TRIM(tide_cpt(jtide)) == 'MU2' ) indx(jtide) = 7 122 if ( TRIM(tide_cpt(jtide)) == 'N2' ) indx(jtide) = 8 123 if ( TRIM(tide_cpt(jtide)) == 'NU2' ) indx(jtide) = 9 124 if ( TRIM(tide_cpt(jtide)) == 'M2' ) indx(jtide) = 10 125 if ( TRIM(tide_cpt(jtide)) == 'L2' ) indx(jtide) = 11 126 if ( TRIM(tide_cpt(jtide)) == 'T2' ) indx(jtide) = 12 127 if ( TRIM(tide_cpt(jtide)) == 'S2' ) indx(jtide) = 13 128 if ( TRIM(tide_cpt(jtide)) == 'K2' ) indx(jtide) = 14 129 if ( TRIM(tide_cpt(jtide)) == 'M4' ) indx(jtide) = 15 130 if (indx(jtide) == 0 ) then 131 if(lwp) write(numout,*) 'WARNING: constitunent', jtide,':', tide_cpt(jtide) & 132 , 'not in standard list' 133 endif 134 enddo 135 ! 136 if ( ntide < 1 ) then 137 if(lwp) write(numout,*) 138 if(lwp) write(numout,*) ' E R R O R : Did not find any tidal components in namelist.' 139 if(lwp) write(numout,*) ' ========== ' 140 if(lwp) write(numout,*) 141 nstop = nstop + 1 142 else 143 if(lwp) write(numout,*) 144 if(lwp) write(numout,*) ntide,' tidal components specified : ' 145 if(lwp) write(numout,*) tide_cpt(1:ntide) 146 if(lwp) write(numout,*) ntide,' phase speeds (deg/hr) : ' 147 if(lwp) write(numout,*) tide_speed(1:ntide) 148 endif 88 ! ! Count number of components specified 89 ntide = jptides_max 90 itide = 1 91 DO WHILE( tide_cpt(itide) /= '' ) 92 ntide = itide 93 itide = itide + 1 94 END DO 95 ! ! find constituents in standard list 96 DO itide = 1, ntide 97 nindx(itide) = 0 98 IF( TRIM( tide_cpt(itide) ) == 'Q1' ) nindx(itide) = 1 99 IF( TRIM( tide_cpt(itide) ) == 'O1' ) nindx(itide) = 2 100 IF( TRIM( tide_cpt(itide) ) == 'P1' ) nindx(itide) = 3 101 IF( TRIM( tide_cpt(itide) ) == 'S1' ) nindx(itide) = 4 102 IF( TRIM( tide_cpt(itide) ) == 'K1' ) nindx(itide) = 5 103 IF( TRIM( tide_cpt(itide) ) == '2N2' ) nindx(itide) = 6 104 IF( TRIM( tide_cpt(itide) ) == 'MU2' ) nindx(itide) = 7 105 IF( TRIM( tide_cpt(itide) ) == 'N2' ) nindx(itide) = 8 106 IF( TRIM( tide_cpt(itide) ) == 'NU2' ) nindx(itide) = 9 107 IF( TRIM( tide_cpt(itide) ) == 'M2' ) nindx(itide) = 10 108 IF( TRIM( tide_cpt(itide) ) == 'L2' ) nindx(itide) = 11 109 IF( TRIM( tide_cpt(itide) ) == 'T2' ) nindx(itide) = 12 110 IF( TRIM( tide_cpt(itide) ) == 'S2' ) nindx(itide) = 13 111 IF( TRIM( tide_cpt(itide) ) == 'K2' ) nindx(itide) = 14 112 IF( TRIM( tide_cpt(itide) ) == 'M4' ) nindx(itide) = 15 113 IF( nindx(itide) == 0 .AND. lwp ) THEN 114 WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list' 115 CALL ctl_warn( ctmp1 ) 116 ENDIF 117 END DO 118 ! ! Parameter control and print 119 IF( ntide < 1 ) THEN 120 CALL ctl_stop( ' Did not find any tidal components in namelist namtide' ) 121 ELSE 122 IF(lwp) WRITE(numout,*) ' Namelist namtide : tidal harmonic forcing at open boundaries' 123 IF(lwp) WRITE(numout,*) ' tidal components specified ', ntide 124 IF(lwp) WRITE(numout,*) ' ', tide_cpt(1:ntide) 125 IF(lwp) WRITE(numout,*) ' associated phase speeds (deg/hr) : ' 126 IF(lwp) WRITE(numout,*) ' ', tide_speed(1:ntide) 127 ENDIF 149 128 150 129 ! Initialisation of tidal harmonics arrays 151 152 130 sshtide(:) = 0.e0 153 utide (:) = 0.e0154 vtide (:) = 0.e0155 131 utide (:) = 0.e0 132 vtide (:) = 0.e0 133 ! 156 134 END SUBROUTINE tide_init 157 135 136 158 137 SUBROUTINE tide_data 138 !!---------------------------------------------------------------------- 139 !! *** SUBROUTINE tide_data *** 140 !! 141 !! ** Purpose : - Read in tidal harmonics data and adjust for the start 142 !! time of the model run. 143 !! 144 !!---------------------------------------------------------------------- 145 INTEGER :: itide, igrd, ib ! dummy loop indices 146 CHARACTER(len=80) :: clfile ! full file name for tidal input file 147 INTEGER :: ipi, ipj, inum, idvar ! temporary integers (netcdf read) 148 INTEGER, DIMENSION(3) :: lendta ! length of data in the file (note may be different from nblendta!) 149 REAL(wp) :: z_arg, z_atde, z_btde, z1t, z2t 150 REAL(wp), DIMENSION(jpbdta,1) :: zdta ! temporary array for data fields 151 REAL(wp), DIMENSION(jptides_max) :: z_vplu, z_ftc 159 152 !!------------------------------------------------------------------------------ 160 !! SUBROUTINE tide_data 161 !! *********************** 162 !! ** Purpose : - Read in tidal harmonics data and adjust for the start time of 163 !! the model run. 164 !! 165 !! History : 166 !! NEMO v2.0 ! 07-01 (D.Storkey) Original 167 !! ! 08-01 (J.Holt) Add date correction. 168 !!------------------------------------------------------------------------------ 169 !! * Local declarations 170 CHARACTER(len=80) :: tidefile ! full file name for tidal input file 171 INTEGER :: jtide, jgrd, jb ! dummy loop indices 172 INTEGER :: ipi, ipj, inum, idvar ! temporary integers (netcdf read) 173 INTEGER, DIMENSION(3) :: lendta ! length of data in the file (note may be 174 ! different from nblendta!) 175 INTEGER :: iday,imonth,iyear 176 REAL(wp) :: arg, atde, btde,z1t,z2t 177 REAL(wp), DIMENSION(jpbdta,1) :: & 178 pdta ! temporary array for data fields 179 REAL(WP), DIMENSION(ntide_max) :: vplu, ftc 180 181 !!------------------------------------------------------------------------------ 182 !! OPA 9.0, LODYC-IPSL (2007) 183 !!------------------------------------------------------------------------------ 184 185 ! 1. Open files and read in tidal forcing data 186 ! -------------------------------------------- 153 154 ! Open files and read in tidal forcing data 155 ! ----------------------------------------- 187 156 188 157 ipj = 1 189 158 190 do jtide = 1 ,ntide 191 192 ! SSH fields 193 ! ---------- 194 195 tidefile = TRIM(filtide)//TRIM(tide_cpt(jtide))//'_grid_T.nc' 196 if(lwp) write(numout,*) 'Reading data from file ',tidefile 197 CALL iom_open( tidefile, inum ) 198 199 jgrd = 1 200 IF ( nblendta(jgrd) .le. 0 ) THEN 201 idvar = iom_varid( inum,'z1' ) 202 if(lwp) write(numout,*) 'iom_file(1)%ndims(idvar) : ',iom_file%ndims(idvar) 203 nblendta(jgrd) = iom_file(inum)%dimsz(1,idvar) 204 WRITE(numout,*) 'Dim size for z1 is ',nblendta(jgrd) 159 DO itide = 1, ntide 160 ! ! SSH fields 161 clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc' 162 IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 163 CALL iom_open( clfile, inum ) 164 igrd = 1 165 IF( nblendta(igrd) <= 0 ) THEN 166 idvar = iom_varid( inum,'z1' ) 167 IF(lwp) WRITE(numout,*) 'iom_file(1)%ndims(idvar) : ',iom_file%ndims(idvar) 168 nblendta(igrd) = iom_file(inum)%dimsz(1,idvar) 169 WRITE(numout,*) 'Dim size for z1 is ', nblendta(igrd) 205 170 ENDIF 206 ipi=nblendta(jgrd) 207 208 CALL iom_get ( inum, jpdom_unknown, 'z1', pdta(1:ipi,1:ipj) ) 209 DO jb=1, nblenrim(jgrd) 210 ssh1(jb,jtide) = pdta(nbmap(jb,jgrd),1) 211 END DO 212 213 CALL iom_get ( inum, jpdom_unknown, 'z2', pdta(1:ipi,1:ipj) ) 214 DO jb=1, nblenrim(jgrd) 215 ssh2(jb,jtide) = pdta(nbmap(jb,jgrd),1) 216 END DO 217 171 ipi = nblendta(igrd) 172 CALL iom_get( inum, jpdom_unknown, 'z1', zdta(1:ipi,1:ipj) ) 173 DO ib = 1, nblenrim(igrd) 174 ssh1(ib,itide) = zdta(nbmap(ib,igrd),1) 175 END DO 176 CALL iom_get( inum, jpdom_unknown, 'z2', zdta(1:ipi,1:ipj) ) 177 DO ib = 1, nblenrim(igrd) 178 ssh2(ib,itide) = zdta(nbmap(ib,igrd),1) 179 END DO 218 180 CALL iom_close( inum ) 219 220 ! U fields 221 ! -------- 222 223 tidefile = TRIM(filtide)//TRIM(tide_cpt(jtide))//'_grid_U.nc' 224 if(lwp) write(numout,*) 'Reading data from file ',tidefile 225 CALL iom_open( tidefile, inum ) 226 227 jgrd = 2 228 IF ( lendta(jgrd) .le. 0 ) THEN 229 idvar = iom_varid( inum,'u1' ) 230 lendta(jgrd) = iom_file(inum)%dimsz(1,idvar) 231 WRITE(numout,*) 'Dim size for u1 is ',lendta(jgrd) 181 ! 182 ! ! U fields 183 clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc' 184 IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 185 CALL iom_open( clfile, inum ) 186 igrd = 2 187 IF( lendta(igrd) <= 0 ) THEN 188 idvar = iom_varid( inum,'u1' ) 189 lendta(igrd) = iom_file(inum)%dimsz(1,idvar) 190 WRITE(numout,*) 'Dim size for u1 is ',lendta(igrd) 232 191 ENDIF 233 ipi=lendta(jgrd) 234 235 CALL iom_get ( inum, jpdom_unknown, 'u1', pdta(1:ipi,1:ipj) ) 236 DO jb=1, nblenrim(jgrd) 237 u1(jb,jtide) = pdta(nbmap(jb,jgrd),1) 238 END DO 239 240 CALL iom_get ( inum, jpdom_unknown, 'u2', pdta(1:ipi,1:ipj) ) 241 DO jb=1, nblenrim(jgrd) 242 u2(jb,jtide) = pdta(nbmap(jb,jgrd),1) 243 END DO 244 192 ipi = lendta(igrd) 193 CALL iom_get( inum, jpdom_unknown, 'u1', zdta(1:ipi,1:ipj) ) 194 DO ib = 1, nblenrim(igrd) 195 u1(ib,itide) = zdta(nbmap(ib,igrd),1) 196 END DO 197 CALL iom_get( inum, jpdom_unknown, 'u2', zdta(1:ipi,1:ipj) ) 198 DO ib = 1, nblenrim(igrd) 199 u2(ib,itide) = zdta(nbmap(ib,igrd),1) 200 END DO 245 201 CALL iom_close( inum ) 246 247 ! V fields 248 ! -------- 249 250 tidefile = TRIM(filtide)//TRIM(tide_cpt(jtide))//'_grid_V.nc' 251 if(lwp) write(numout,*) 'Reading data from file ',tidefile 252 CALL iom_open( tidefile, inum ) 253 254 jgrd = 3 255 IF ( lendta(jgrd) .le. 0 ) THEN 256 idvar = iom_varid( inum,'v1' ) 257 lendta(jgrd) = iom_file(inum)%dimsz(1,idvar) 258 WRITE(numout,*) 'Dim size for v1 is ',lendta(jgrd) 202 ! 203 ! ! V fields 204 clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc' 205 if(lwp) write(numout,*) 'Reading data from file ', clfile 206 CALL iom_open( clfile, inum ) 207 igrd = 3 208 IF( lendta(igrd) <= 0 ) THEN 209 idvar = iom_varid( inum,'v1' ) 210 lendta(igrd) = iom_file(inum)%dimsz(1,idvar) 211 WRITE(numout,*) 'Dim size for v1 is ', lendta(igrd) 259 212 ENDIF 260 ipi=lendta(jgrd) 261 262 CALL iom_get ( inum, jpdom_unknown, 'v1', pdta(1:ipi,1:ipj) ) 263 DO jb=1, nblenrim(jgrd) 264 v1(jb,jtide) = pdta(nbmap(jb,jgrd),1) 265 END DO 266 267 CALL iom_get ( inum, jpdom_unknown, 'v2', pdta(1:ipi,1:ipj) ) 268 DO jb=1, nblenrim(jgrd) 269 v2(jb,jtide) = pdta(nbmap(jb,jgrd),1) 270 END DO 271 213 ipi = lendta(igrd) 214 CALL iom_get( inum, jpdom_unknown, 'v1', zdta(1:ipi,1:ipj) ) 215 DO ib = 1, nblenrim(igrd) 216 v1(ib,itide) = zdta(nbmap(ib,igrd),1) 217 END DO 218 CALL iom_get( inum, jpdom_unknown, 'v2', zdta(1:ipi,1:ipj) ) 219 DO ib=1, nblenrim(igrd) 220 v2(ib,itide) = zdta(nbmap(ib,igrd),1) 221 END DO 272 222 CALL iom_close( inum ) 273 enddo ! end loop on tidal components 274 ! 275 ! correct for date factors 276 ! 277 if( ln_tide_date ) then 278 279 280 ! Calculate date corrects for 15 standard consituents 281 iyear = int(ndate0 / 10000 ) ! initial year 282 imonth = int((ndate0 - iyear * 10000 ) / 100 ) ! initial month 283 iday = int(ndate0 - iyear * 10000 - imonth * 100) ! initial day betwen 1 and 30 284 285 call uvset(0,iday,imonth,iyear,ftc,vplu) 286 ! 287 if(lwp) write(numout,*) 'Correcting tide for date:',iday,imonth,iyear 288 289 do jtide = 1, ntide 290 ! 291 if(indx(jtide) .ne. 0) then 292 arg=3.14159265d0*vplu(indx(jtide))/180.0d0 293 atde=ftc(indx(jtide))*cos(arg) 294 btde=ftc(indx(jtide))*sin(arg) 295 if(lwp) then 296 write(numout,'(2i5,8f10.6)') jtide,indx(jtide),tide_speed(jtide), & 297 ftc(indx(jtide)),vplu(indx(jtide)) 298 endif 299 else 300 atde = 1.0_wp 301 btde = 0.0_wp 302 endif 303 ! elevation 304 jgrd = 1 305 do jb=1, nblenrim(jgrd) 306 z1t=atde*ssh1(jb,jtide)+btde*ssh2(jb,jtide) 307 z2t=atde*ssh2(jb,jtide)-btde*ssh1(jb,jtide) 308 ssh1(jb,jtide) = z1t 309 ssh2(jb,jtide) = z2t 310 end do 311 ! u 312 jgrd = 2 313 do jb=1, nblenrim(jgrd) 314 z1t=atde*u1(jb,jtide)+btde*u2(jb,jtide) 315 z2t=atde*u2(jb,jtide)-btde*u1(jb,jtide) 316 u1(jb,jtide) = z1t 317 u2(jb,jtide) = z2t 318 end do 319 ! v 320 jgrd = 3 321 do jb=1, nblenrim(jgrd) 322 z1t=atde*v1(jb,jtide)+btde*v2(jb,jtide) 323 z2t=atde*v2(jb,jtide)-btde*v1(jb,jtide) 324 v1(jb,jtide) = z1t 325 v2(jb,jtide) = z2t 326 end do 327 328 enddo ! end loop on tidal components 329 endif ! date correction 330 223 ! 224 END DO ! end loop on tidal components 225 226 IF( ln_tide_date ) THEN ! correct for date factors 227 228 !! used nmonth, nyear and nday from daymod.... 229 ! Calculate date corrects for 15 standard consituents 230 ! This is the initialisation step, so nday, nmonth, nyear are the 231 ! initial date/time of the integration. 232 233 CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu ) 234 235 IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear 236 237 DO itide = 1, ntide ! loop on tidal components 238 ! 239 IF( nindx(itide) /= 0 ) THEN 240 !!gm use rpi and rad global variable 241 z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0 242 z_atde=z_ftc(nindx(itide))*cos(z_arg) 243 z_btde=z_ftc(nindx(itide))*sin(z_arg) 244 IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), tide_speed(itide), & 245 & z_ftc(nindx(itide)), z_vplu(nindx(itide)) 246 ELSE 247 z_atde = 1.0_wp 248 z_btde = 0.0_wp 249 ENDIF 250 ! ! elevation 251 igrd = 1 252 DO ib = 1, nblenrim(igrd) 253 z1t = z_atde * ssh1(ib,itide) + z_btde * ssh2(ib,itide) 254 z2t = z_atde * ssh2(ib,itide) - z_btde * ssh1(ib,itide) 255 ssh1(ib,itide) = z1t 256 ssh2(ib,itide) = z2t 257 END DO 258 ! ! u 259 igrd = 2 260 DO ib = 1, nblenrim(igrd) 261 z1t = z_atde * u1(ib,itide) + z_btde * u2(ib,itide) 262 z2t = z_atde * u2(ib,itide) - z_btde * u1(ib,itide) 263 u1(ib,itide) = z1t 264 u2(ib,itide) = z2t 265 END DO 266 ! ! v 267 igrd = 3 268 DO ib = 1, nblenrim(igrd) 269 z1t = z_atde * v1(ib,itide) + z_btde * v2(ib,itide) 270 z2t = z_atde * v2(ib,itide) - z_btde * v1(ib,itide) 271 v1(ib,itide) = z1t 272 v2(ib,itide) = z2t 273 END DO 274 ! 275 END DO ! end loop on tidal components 276 ! 277 ENDIF ! date correction 278 ! 331 279 END SUBROUTINE tide_data 332 280 281 333 282 SUBROUTINE tide_update ( kt, jit ) 334 !!---------------------------------------------------------------------- --------335 !! SUBROUTINE tide_update336 !! ************************283 !!---------------------------------------------------------------------- 284 !! *** SUBROUTINE tide_update *** 285 !! 337 286 !! ** Purpose : - Add tidal forcing to sshbdy, ubtbdy and vbtbdy arrays. 338 287 !! 339 !! 340 !! History : 341 !! NEMO v2.0 ! 06-12 (D.Storkey) Original 342 !!------------------------------------------------------------------------------ 343 !! * Arguments 288 !!---------------------------------------------------------------------- 344 289 INTEGER, INTENT( in ) :: kt ! Main timestep counter 290 !!gm doctor jit ==> kit 345 291 INTEGER, INTENT( in ) :: jit ! Barotropic timestep counter (for timesplitting option) 346 347 !! * Local declarations 348 INTEGER :: jtide, jgrd, jb ! dummy loop indices 349 REAL(wp) :: arg, sarg 350 REAL(wp), DIMENSION(ntide_max) :: sist, cost 351 352 !!------------------------------------------------------------------------------ 353 !! OPA 9.0, LODYC-IPSL (2003) 354 !!------------------------------------------------------------------------------ 292 !! 293 INTEGER :: itide, igrd, ib ! dummy loop indices 294 REAL(wp) :: z_arg, z_sarg ! 295 REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost 296 !!---------------------------------------------------------------------- 355 297 356 298 ! Note tide phase speeds are in deg/hour, so we need to convert the 357 299 ! elapsed time in seconds to hours by dividing by 3600.0 358 if ( jit .eq. 0 ) then 359 arg = kt*rdt*rad/3600.0 360 else ! we are in a barotropic subcycle (for timesplitting option) 361 arg = ( (kt-1)*rdt + jit*rdtbt ) * rad/3600.0 362 endif 363 364 do jtide = 1,ntide 365 sarg = arg*tide_speed(jtide) 366 cost(jtide) = cos(sarg) 367 sist(jtide) = sin(sarg) 368 enddo 369 370 !! summing of tidal constituents into BDY arrays 371 300 IF( jit == 0 ) THEN 301 z_arg = kt * rdt * rad / 3600.0 302 ELSE ! we are in a barotropic subcycle (for timesplitting option) 303 z_arg = ( (kt-1) * rdt + jit * rdtbt ) * rad / 3600.0 304 ENDIF 305 306 DO itide = 1, ntide 307 z_sarg = z_arg * tide_speed(itide) 308 z_cost(itide) = COS( z_sarg ) 309 z_sist(itide) = SIN( z_sarg ) 310 END DO 311 312 ! summing of tidal constituents into BDY arrays 372 313 sshtide(:) = 0.0 373 utide(:) = 0.0 374 vtide(:) = 0.0 375 376 do jtide = 1 ,ntide 377 jgrd=1 !: SSH on tracer grid. 378 do jb = 1, nblenrim(jgrd) 379 sshtide(jb) =sshtide(jb)+ ssh1(jb,jtide)*cost(jtide) + ssh2(jb,jtide)*sist(jtide) 380 ! if(lwp) write(numout,*) 'z',jb,jtide,sshtide(jb), ssh1(jb,jtide),ssh2(jb,jtide) 381 enddo 382 jgrd=2 !: U grid 383 do jb=1, nblenrim(jgrd) 384 utide(jb) = utide(jb)+ u1(jb,jtide)*cost(jtide) + u2(jb,jtide)*sist(jtide) 385 ! if(lwp) write(numout,*) 'u',jb,jtide,utide(jb), u1(jb,jtide),u2(jb,jtide) 386 387 enddo 388 jgrd=3 !: V grid 389 do jb=1, nblenrim(jgrd) 390 vtide(jb) = vtide(jb)+ v1(jb,jtide)*cost(jtide) + v2(jb,jtide)*sist(jtide) 391 ! if(lwp) write(numout,*) 'v',jb,jtide,vtide(jb), v1(jb,jtide),v2(jb,jtide) 392 393 enddo 394 enddo 395 314 utide (:) = 0.0 315 vtide (:) = 0.0 316 ! 317 DO itide = 1, ntide 318 igrd=1 ! SSH on tracer grid. 319 DO ib = 1, nblenrim(igrd) 320 sshtide(ib) =sshtide(ib)+ ssh1(ib,itide)*z_cost(itide) + ssh2(ib,itide)*z_sist(itide) 321 ! if(lwp) write(numout,*) 'z',ib,itide,sshtide(ib), ssh1(ib,itide),ssh2(ib,itide) 322 END DO 323 igrd=2 ! U grid 324 DO ib=1, nblenrim(igrd) 325 utide(ib) = utide(ib)+ u1(ib,itide)*z_cost(itide) + u2(ib,itide)*z_sist(itide) 326 ! if(lwp) write(numout,*) 'u',ib,itide,utide(ib), u1(ib,itide),u2(ib,itide) 327 END DO 328 igrd=3 ! V grid 329 DO ib=1, nblenrim(igrd) 330 vtide(ib) = vtide(ib)+ v1(ib,itide)*z_cost(itide) + v2(ib,itide)*z_sist(itide) 331 ! if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), v1(ib,itide),v2(ib,itide) 332 END DO 333 END DO 334 ! 396 335 END SUBROUTINE tide_update 397 336 398 ! 399 ! 400 SUBROUTINE uvset (ihs,iday,imnth,iyr,f,vplu)401 !!---------------------------------------------------------------------- --------402 !! SUBROUTINE uvset403 !! ************************337 !!gm doctor naming of dummy argument variables!!! and all local variables 338 339 SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu ) 340 !!---------------------------------------------------------------------- 341 !! *** SUBROUTINE uvset *** 342 !! 404 343 !! ** Purpose : - adjust tidal forcing for date factors 405 344 !! 406 !! 407 !! History : 408 !! 409 !! Origins POLCOMS v6.3 2007 410 !! NEMO v2.3 ! Jason Holt 411 !!------------------------------------------------------------------------------ 345 !!---------------------------------------------------------------------- 412 346 implicit none 413 !! * Arguments414 347 INTEGER, INTENT( in ) :: ihs ! Start time hours 415 348 INTEGER, INTENT( in ) :: iday ! start time days 416 INTEGER, INTENT( in ) :: imnth ! start time month 417 INTEGER, INTENT( in ) :: iyr ! start time year 418 INTEGER, PARAMETER :: nc =15 ! maximum number of constituents 419 420 REAL(WP) :: f(nc) ! nodal correction 421 REAL(WP) :: vplu(nc) ! phase correction 422 ! 423 INTEGER :: year,vd,ivdy,ndc,i,k 424 REAL(WP) :: u(nc),v(nc),zig(nc),rtd 425 REAL(WP) :: ss,h,p,en,p1 349 INTEGER, INTENT( in ) :: imnth ! start time month 350 INTEGER, INTENT( in ) :: iyr ! start time year 351 !! 352 !!gm nc is jptides_max ???? 353 INTEGER , PARAMETER :: nc =15 ! maximum number of constituents 426 354 CHARACTER(len=8), DIMENSION(nc) :: cname 427 ! 428 !!------------------------------------------------------------------------------ 429 !! OPA 9.0, LODYC-IPSL (2007) 430 !!------------------------------------------------------------------------------ 431 432 data cname/ 'q1', 'o1', 'p1', 's1', 'k1', & 433 '2n2','mu2', 'n2','nu2', 'm2', 'l2', 't2', 's2', 'k2', & 434 'm4'/ 435 data zig/.2338507481,.2433518789,.2610826055,.2617993878 & 436 , .2625161701 & 437 , .4868657873,.4881373225,.4963669182,.4976384533 & 438 , .5058680490,.5153691799,.5228820265,.5235987756 & 439 , .5250323419 & 440 , 1.011736098/ 355 INTEGER :: year, vd, ivdy, ndc, i, k 356 REAL(wp) :: ss, h, p, en, p1, rtd 357 REAL(wp), DIMENSION(nc) :: f ! nodal correction 358 REAL(wp), DIMENSION(nc) :: z_vplu ! phase correction 359 REAL(wp), DIMENSION(nc) :: u, v, zig 360 !! 361 DATA cname/ 'q1' , 'o1' , 'p1' , 's1' , 'k1' , & 362 & '2n2' , 'mu2' , 'n2' , 'nu2' , 'm2' , & 363 & 'l2' , 't2' , 's2' , 'k2' , 'm4' / 364 DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878, .2625161701, & 365 & .4868657873, .4881373225, .4963669182, .4976384533, .5058680490, & 366 & .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098 / 367 !!---------------------------------------------------------------------- 441 368 ! 442 369 ! ihs - start time gmt on ... 443 370 ! iday/imnth/iyr - date e.g. 12/10/87 444 371 ! 445 callvday(iday,imnth,iyr,ivdy)446 vd =ivdy372 CALL vday(iday,imnth,iyr,ivdy) 373 vd = ivdy 447 374 ! 448 375 !rp note change of year number for d. blackman shpen … … 453 380 !.....year = year of required data 454 381 !.....vd = day of required data..set up for 0000gmt day year 455 !456 382 ndc = nc 457 383 !.....ndc = number of constituents allowed 458 ! 459 rtd = 360.0/6.2831852 460 do i = 1,ndc 461 zig(i) = zig(i)*rtd 462 ! sigo(i)= zig(i) 463 enddo 464 ! 465 if(year == 0) then 466 return 467 endif 384 !!gm use rpi ? 385 rtd = 360.0 / 6.2831852 386 DO i = 1, ndc 387 zig(i) = zig(i)*rtd 388 ! sigo(i)= zig(i) 389 END DO 390 391 !!gm try to avoid RETURN in F90 392 IF( year == 0 ) RETURN 468 393 469 call shpen (year,vd,ss,h,p,en,p1) 470 call ufset (p,en,u,f) 471 call vset (ss,h,p,en,p1,v) 472 ! 473 do k=1,nc 474 475 vplu(k)=v(k)+u(k) 476 vplu(k)=vplu(k)+dble(ihs)*zig(k) 477 ! 478 do while ( vplu(k) < 0 ) 479 vplu(k) = vplu(k) + 360.0 480 enddo 481 ! 482 do while (vplu(k) > 360. ) 483 vplu(k) = vplu(k) - 360.0 484 enddo 485 ! 486 enddo 487 ! 394 CALL shpen( year, vd, ss, h , p , en, p1 ) 395 CALL ufset( p , en, u , f ) 396 CALL vset ( ss , h , p , en, p1, v ) 397 ! 398 DO k = 1, nc 399 z_vplu(k) = v(k) + u(k) 400 z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k) 401 DO WHILE( z_vplu(k) < 0 ) 402 z_vplu(k) = z_vplu(k) + 360.0 403 END DO 404 DO WHILE( z_vplu(k) > 360. ) 405 z_vplu(k) = z_vplu(k) - 360.0 406 END DO 407 END DO 408 ! 488 409 END SUBROUTINE uvset 489 410 490 411 491 SUBROUTINE vday( iday,imnth,iy,ivdy)492 !!---------------------------------------------------------------------- --------493 !! SUBROUTINE vday494 !! ****************412 SUBROUTINE vday( iday, imnth, iy, ivdy ) 413 !!---------------------------------------------------------------------- 414 !! *** SUBROUTINE vday *** 415 !! 495 416 !! ** Purpose : - adjust tidal forcing for date factors 496 417 !! 497 !! 498 !! History :499 !!500 !! Origins POLCOMS v6.3 2007501 !! NEMO v2.3 ! Jason Holt418 !!---------------------------------------------------------------------- 419 INTEGER, INTENT(in ) :: iday, imnth, iy ! ???? 420 INTEGER, INTENT( out) :: ivdy ! ??? 421 !! 422 INTEGER :: iyr 502 423 !!------------------------------------------------------------------------------ 503 implicit none 504 ! 505 ! subroutine arguments 506 integer :: iday,imnth,iy,ivdy 507 ! 508 ! local variables 509 integer iyr 510 !!------------------------------------------------------------------------------ 511 !! NEMO 2.3, LODYC-IPSL (2008) 512 !!------------------------------------------------------------------------------ 513 514 ! ================================================================= 515 ! 516 ! calculate day number in year from day/month/year 517 ! 518 ! ================================================================= 424 425 !!gm nday_year in day mode is the variable compiuted here, no? 426 !!gm nday_year , & !: curent day counted from jan 1st of the current year 427 428 !calculate day number in year from day/month/year 519 429 if(imnth.eq.1) ivdy=iday 520 430 if(imnth.eq.2) ivdy=iday+31 … … 529 439 if(imnth.eq.11) ivdy=iday+304 530 440 if(imnth.eq.12) ivdy=iday+334 531 iyr=iy441 iyr=iy 532 442 if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1 533 443 if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1 534 444 if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1 535 536 RETURN 537 END SUBROUTINE vday 538 539 SUBROUTINE shpen (year,vd,s,h,p,en,p1) 540 !!------------------------------------------------------------------------------ 541 !! SUBROUTINE shpen 542 !! ***************** 445 ! 446 END SUBROUTINE vday 447 448 !!doctor norme for dummy arguments 449 450 SUBROUTINE shpen( year, vd, s, h, p, en, p1 ) 451 !!---------------------------------------------------------------------- 452 !! *** SUBROUTINE shpen *** 453 !! 543 454 !! ** Purpose : - calculate astronomical arguments for tides 544 455 !! this version from d. blackman 30 nove 1990 545 !! 546 !! History : 547 !! 548 !! Origins POLCOMS v6.3 2007 549 !! NEMO v2.3 ! Jason Holt 550 !!------------------------------------------------------------------------------ 551 implicit none 552 553 ! subroutine arguments 554 integer ::year,vd 555 real(wp) :: s,h,p,en,p1 556 557 ! local variables 558 559 integer yr,ilc,icent,it,iday,ild,ipos,nn,iyd 456 !! 457 !!---------------------------------------------------------------------- 458 !!gm add INTENT in, out or inout.... DOCTOR name.... 459 !!gm please do not use variable name with a single letter: impossible to search in a code 460 INTEGER :: year,vd 461 REAL(wp) :: s,h,p,en,p1 462 !! 463 INTEGER :: yr,ilc,icent,it,iday,ild,ipos,nn,iyd 560 464 REAL(wp) :: cycle,t,td,delt(84),delta,deltat 561 562 !!------------------------------------------------------------------------------ 563 !! NEMO 2.3, LODYC-IPSL (2008) 564 !!------------------------------------------------------------------------------ 565 566 data delt /-5.04,-3.90,-2.87,-0.58,0.71,1.80, & 567 3.08, 4.63, 5.86, 7.21, 8.58,10.50,12.10, & 568 12.49,14.41,15.59,15.81,17.52,19.01,18.39, & 569 19.55,20.36,21.01,21.81,21.76,22.35,22.68, & 570 22.94,22.93,22.69,22.94,23.20,23.31,23.63, & 571 23.47,23.68,23.62,23.53,23.59,23.99,23.80, & 572 24.20,24.99,24.97,25.72,26.21,26.37,26.89, & 573 27.68,28.13,28.94,29.42,29.66,30.29,30.96, & 574 31.09,31.59,31.52,31.92,32.45,32.91,33.39, & 575 33.80,34.23,34.73,35.40,36.14,36.99,37.87, & 576 38.75,39.70,40.70,41.68,42.82,43.96,45.00, & 577 45.98,47.00,48.03,49.10,50.10,50.97,51.81, & 578 52.57/ 579 ! 465 !! 466 DATA delt /-5.04, -3.90, -2.87, -0.58, 0.71, 1.80, & 467 & 3.08, 4.63, 5.86, 7.21, 8.58, 10.50, 12.10, & 468 & 12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39, & 469 & 19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68, & 470 & 22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63, & 471 & 23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80, & 472 & 24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89, & 473 & 27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96, & 474 & 31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39, & 475 & 33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87, & 476 & 38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00, & 477 & 45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81, & 478 & 52.57 / 479 !!---------------------------------------------------------------------- 480 580 481 cycle = 360.0 581 482 ilc = 0 582 icent = year /100583 yr = year - icent *100483 icent = year / 100 484 yr = year - icent * 100 584 485 t = icent - 20 585 486 ! … … 588 489 ! see notes by cartwright 589 490 ! 491 !!gm old coding style, use CASE instead and avoid GOTO (obsolescence in fortran 90) 492 !!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 590 493 it = icent - 20 591 494 if (it) 1,2,2 … … 602 505 td = 0.0 603 506 ! 507 !!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 604 508 if (yr) 4,5,4 605 509 ! … … 623 527 7 deltat = delta * 1.0e-6 624 528 ! 529 !!gm precision of the computation : example for s it should be replace by: 530 !!gm s = 218.3165 + (481267.8813 - 0.0016*t)*t + 152.0*deltat ==> more precise modify the last digits results 625 531 s = 218.3165 + 481267.8813*t - 0.0016*t*t + 152.0*deltat 626 h = 280.4661 + 36000.7698 *t + 0.0003*t*t +11.0*deltat627 p = 83.3535 + 4069.0139 *t - 0.0103*t*t +deltat628 en = 234.9555 + 1934.1363 *t - 0.0021*t*t +deltat629 p1 = 282.9384 + 1.7195 *t + 0.0005*t*t630 !631 nn = s /cycle632 s = s - nn *cycle633 if ( s .lt. 0.0) s = s+cycle634 !635 nn = h /cycle636 h = h -cycle*nn637 if (h .lt. 0.0) h = h+cycle638 !639 nn = p /cycle640 p = p - cycle*nn641 if (p .lt. 0.0) p = p+cycle642 !643 nn = en /cycle644 en = en -cycle*nn645 if(en .lt. 0.0)en = en + cycle532 h = 280.4661 + 36000.7698 *t + 0.0003*t*t + 11.0*deltat 533 p = 83.3535 + 4069.0139 *t - 0.0103*t*t + deltat 534 en = 234.9555 + 1934.1363 *t - 0.0021*t*t + deltat 535 p1 = 282.9384 + 1.7195 *t + 0.0005*t*t 536 ! 537 nn = s / cycle 538 s = s - nn * cycle 539 IF( s < 0.e0 ) s = s + cycle 540 ! 541 nn = h / cycle 542 h = h - cycle * nn 543 IF( h < 0.e0 ) h = h + cycle 544 ! 545 nn = p / cycle 546 p = p - cycle * nn 547 IF( p < 0.e0) p = p + cycle 548 ! 549 nn = en / cycle 550 en = en - cycle * nn 551 IF( en < 0.e0 ) en = en + cycle 646 552 en = cycle - en 647 ! 648 nn = p1/cycle 649 p1 = p1 - nn*cycle 650 ! 651 RETURN 652 ! 653 END SUBROUTINE shpen 654 655 656 SUBROUTINE ufset (p,cn,b,a) 657 !!------------------------------------------------------------------------------ 658 !! SUBROUTINE ufset 659 !! ***************** 553 ! 554 nn = p1 / cycle 555 p1 = p1 - nn * cycle 556 ! 557 END SUBROUTINE shpen 558 559 560 SUBROUTINE ufset( p, cn, b, a ) 561 !!---------------------------------------------------------------------- 562 !! *** SUBROUTINE ufset *** 563 !! 660 564 !! ** Purpose : - calculate nodal parameters for the tides 661 565 !! 662 !! 663 !! History : 664 !! 665 !! Origins POLCOMS v6.3 2007 666 !! NEMO v2.3 ! Jason Holt 667 !!------------------------------------------------------------------------------ 668 669 implicit none 566 !!---------------------------------------------------------------------- 567 !!gm doctor naming of dummy argument variables!!! and all local variables 568 !!gm nc is jptides_max ???? 670 569 integer nc 671 570 parameter (nc=15) 672 ! subroutine arguments 673 real(wp) p,cn 674 ! 675 ! 676 ! local variables 677 real(wp) :: w1,w2,w3,w4,w5,w6,w7,w8,nw,pw,rad 678 real(wp) :: a(nc),b(nc) 679 integer ndc,k 680 ! 681 !!------------------------------------------------------------------------------ 682 !! NEMO 2.3, LODYC-IPSL (2008) 683 !!------------------------------------------------------------------------------ 684 685 ndc=nc 686 ! 571 REAL(wp) p,cn 572 !! 573 574 !!gm rad is already a public variable defined in phycst.F90 .... ==> doctor norme local real start with "z" 575 REAL(wp) :: w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad 576 REAL(wp) :: a(nc), b(nc) 577 INTEGER :: ndc, k 578 !!---------------------------------------------------------------------- 579 580 ndc = nc 581 687 582 ! a=f , b =u 688 583 ! t is zero as compared to tifa. 584 !! use rad defined in phycst (i.e. add a USE phycst at the begining of the module 689 585 rad = 6.2831852d0/360.0 690 pw = p *rad691 nw = cn *rad692 w1 = cos( nw)693 w2 = cos( 2*nw)694 w3 = cos( 3*nw)695 w4 = sin( nw)696 w5 = sin( 2*nw)697 w6 = sin( 3*nw)698 w7 = 1 -0.2505*cos(2*pw) -0.1102*cos(2*pw-nw)&699 -0.156*cos(2*pw-2*nw) -0.037*cos(nw)700 w8 = -0.2505*sin(2*pw) -0.1102*sin(2*pw-nw)&701 -0.0156*sin(2*pw-2*nw) -0.037*sin(nw)702 !703 a(1) = 1.0089 +0.1871*w1-0.0147*w2+0.0014*w3704 b(1) = 0.1885*w4 - 0.0234*w5+.0033*w6705 ! q1586 pw = p * rad 587 nw = cn * rad 588 w1 = cos( nw ) 589 w2 = cos( 2*nw ) 590 w3 = cos( 3*nw ) 591 w4 = sin( nw ) 592 w5 = sin( 2*nw ) 593 w6 = sin( 3*nw ) 594 w7 = 1. - 0.2505 * COS( 2*pw ) - 0.1102 * COS(2*pw-nw ) & 595 & - 0.156 * COS( 2*pw-2*nw ) - 0.037 * COS( nw ) 596 w8 = - 0.2505 * SIN( 2*pw ) - 0.1102 * SIN(2*pw-nw ) & 597 & - 0.0156 * SIN( 2*pw-2*nw ) - 0.037 * SIN( nw ) 598 ! 599 a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3 600 b(1) = 0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6 601 ! q1 706 602 a(2) = a(1) 707 603 b(2) = b(1) 708 ! o1604 ! o1 709 605 a(3) = 1.0 710 606 b(3) = 0.0 711 ! p1607 ! p1 712 608 a(4) = 1.0 713 609 b(4) = 0.0 714 ! s1610 ! s1 715 611 a(5) = 1.0060+0.1150*w1- 0.0088*w2 +0.0006*w3 716 612 b(5) = -0.1546*w4 + 0.0119*w5 -0.0012*w6 717 ! k1613 ! k1 718 614 a(6) =1.0004 -0.0373*w1+ 0.0002*w2 719 615 b(6) = -0.0374*w4 720 ! 2n2616 ! 2n2 721 617 a(7) = a(6) 722 618 b(7) = b(6) 723 ! mu2619 ! mu2 724 620 a(8) = a(6) 725 621 b(8) = b(6) 726 ! n2622 ! n2 727 623 a(9) = a(6) 728 624 b(9) = b(6) 729 ! nu2625 ! nu2 730 626 a(10) = a(6) 731 627 b(10) = b(6) 732 ! m2 733 a(11) = sqrt(w7*w7+w8*w8) 734 b(11) = atan(w8/w7) 735 if(w7.lt.0) b(11) = b(11) + 3.141992 736 ! l2 628 ! m2 629 a(11) = SQRT( w7 * w7 + w8 * w8 ) 630 b(11) = ATAN( w8 / w7 ) 631 !!gmuse rpi instead of 3.141992 ??? true pi is rpi=3.141592653589793_wp ..... ???? 632 IF( w7 < 0.e0 ) b(11) = b(11) + 3.141992 633 ! l2 737 634 a(12) = 1.0 738 635 b(12) = 0.0 739 ! t2636 ! t2 740 637 a(13)= a(12) 741 638 b(13)= b(12) 742 ! s2639 ! s2 743 640 a(14) = 1.0241+0.2863*w1+0.0083*w2 -0.0015*w3 744 641 b(14) = -0.3096*w4 + 0.0119*w5 - 0.0007*w6 745 ! k2642 ! k2 746 643 a(15) = a(6)*a(6) 747 644 b(15) = 2*b(6) 748 ! m4 749 ! 750 do 40 k = 1,ndc 751 b(k) = b(k)/rad 752 32 if (b(k)) 34,35,35 753 34 b(k) = b(k) + 360.0 754 go to 32 755 35 if (b(k)-360.0) 40,37,37 756 37 b(k) = b(k)-360.0 757 go to 35 645 ! m4 646 !!gm old coding, remove GOTO and label of lines 647 !!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 648 DO 40 k = 1,ndc 649 b(k) = b(k)/rad 650 32 if (b(k)) 34,35,35 651 34 b(k) = b(k) + 360.0 652 go to 32 653 35 if (b(k)-360.0) 40,37,37 654 37 b(k) = b(k)-360.0 655 go to 35 758 656 40 continue 759 RETURN 760 END SUBROUTINE ufset 761 762 SUBROUTINE vset( s,h,p,en,p1,v) 763 !!------------------------------------------------------------------------------ 764 !! SUBROUTINE vset 765 !! **************** 657 ! 658 END SUBROUTINE ufset 659 660 661 SUBROUTINE vset( s,h,p,en,p1,v) 662 !!---------------------------------------------------------------------- 663 !! *** SUBROUTINE vset *** 664 !! 766 665 !! ** Purpose : - calculate tidal phases for 0000gmt on start day of run 767 666 !! 768 !! 769 !! History : 770 !! 771 !! Origins POLCOMS v6.3 2007 772 !! NEMO v2.3 ! Jason Holt 773 !!------------------------------------------------------------------------------ 774 implicit none 667 !!---------------------------------------------------------------------- 668 !!gm doctor naming of dummy argument variables!!! and all local variables 669 !!gm nc is jptides_max ???? 670 !!gm en argument is not used: suppress it ? 775 671 integer nc 776 672 parameter (nc=15) 777 ! subroutine arguments778 673 real(wp) s,h,p,en,p1 779 674 real(wp) v(nc) 780 ! 781 integer ndc,k 782 !!------------------------------------------------------------------------------ 783 !! NEMO 2.3, LODYC-IPSL (2008) 784 !!------------------------------------------------------------------------------ 675 !! 676 integer ndc, k 677 !!---------------------------------------------------------------------- 785 678 786 679 ndc = nc 787 ! v s are computed here. 788 v(1) =-3*s +h +p +270 789 ! q1 790 v(2) =-2*s +h +270.0 791 ! o1 792 v(3) =-h +270 793 ! p1 794 v(4) =180 795 ! s1 796 v(5) =h +90.0 797 ! k1 798 v(6) =-4*s +2*h +2*p 799 ! 2n2 800 v(7) =-4*(s-h) 801 ! mu2 802 v(8) =-3*s +2*h +p 803 ! n2 804 v(9) =-3*s +4*h -p 805 ! mu2 806 v(10) =-2*s +2*h 807 ! m2 808 v(11) =-s +2*h -p +180 809 ! l2 810 v(12) =-h +p1 811 ! t2 812 v(13) =0 813 ! s2 814 v(14) =h+h 815 ! k2 816 v(15) =2*v(10) 817 ! m4 818 ! 680 ! v s are computed here. 681 v(1) =-3*s +h +p +270 ! Q1 682 v(2) =-2*s +h +270.0 ! O1 683 v(3) =-h +270 ! P1 684 v(4) =180 ! S1 685 v(5) =h +90.0 ! K1 686 v(6) =-4*s +2*h +2*p ! 2N2 687 v(7) =-4*(s-h) ! MU2 688 v(8) =-3*s +2*h +p ! N2 689 v(9) =-3*s +4*h -p ! MU2 690 v(10) =-2*s +2*h ! M2 691 v(11) =-s +2*h -p +180 ! L2 692 v(12) =-h +p1 ! T2 693 v(13) =0 ! S2 694 v(14) =h+h ! K2 695 v(15) =2*v(10) ! M4 696 ! 697 !!gm old coding, remove GOTO and label of lines 698 !!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 819 699 do 72 k = 1, ndc 820 69 if (v(k) )70,71,71821 70 v(k)=v(k)+360.0822 go to 69823 71 if ( v(k) - 360.0)72,73,73824 73 v(k)=v(k)-360.0825 go to 71700 69 if( v(k) ) 70,71,71 701 70 v(k) = v(k)+360.0 702 go to 69 703 71 if( v(k) - 360.0 ) 72,73,73 704 73 v(k) = v(k)-360.0 705 go to 71 826 706 72 continue 827 RETURN828 707 ! 708 END SUBROUTINE vset 829 709 830 710 #else 831 !!================================================================================= 832 !! *** MODULE bdytides *** 833 !!================================================================================= 834 835 LOGICAL, PUBLIC, PARAMETER :: lk_bdy_tides = .FALSE. !: tidal forcing at boundaries. 836 837 CHARACTER(len=80), PUBLIC :: & 838 filtide !: Filename root for tidal input files 839 840 CHARACTER(len=4), PUBLIC, DIMENSION(1) :: & 841 tide_cpt !: Names of tidal components used. 711 !!---------------------------------------------------------------------- 712 !! Dummy module NO Unstruct Open Boundary Conditions for tides 713 !!---------------------------------------------------------------------- 714 !!gm are you sure we need to define filtide and tide_cpt ? 715 CHARACTER(len=80), PUBLIC :: filtide !: Filename root for tidal input files 716 CHARACTER(len=4 ), PUBLIC, DIMENSION(1) :: tide_cpt !: Names of tidal components used. 842 717 843 718 CONTAINS 844 845 SUBROUTINE tide_init 846 ! No tidal forcing at boundaries ==> empty routine 719 SUBROUTINE tide_init ! Empty routine 847 720 END SUBROUTINE tide_init 848 849 SUBROUTINE tide_data 850 ! No tidal forcing at boundaries ==> empty routine 721 SUBROUTINE tide_data ! Empty routine 851 722 END SUBROUTINE tide_data 852 853 SUBROUTINE tide_update( kt, jit ) 854 INTEGER :: kt, jit 855 WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, jit 856 857 ! No tidal forcing at boundaries ==> empty routine 723 SUBROUTINE tide_update( kt, kit ) ! Empty routine 724 WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, kit 858 725 END SUBROUTINE tide_update 859 860 726 #endif 861 727 728 !!====================================================================== 862 729 END MODULE bdytides -
trunk/NEMO/OPA_SRC/BDY/bdytra.F90
r911 r1125 1 1 MODULE bdytra 2 !!====================================================================== ===========2 !!====================================================================== 3 3 !! *** MODULE bdytra *** 4 4 !! Ocean tracers: Flow Relaxation Scheme of tracers on each open boundary 5 !!================================================================================= 5 !!====================================================================== 6 !! History : 1.0 ! 2005-01 (J. Chanut, A. Sellar) Original code 7 !! 3.0 ! 2008-04 (NEMO team) add in the reference version 8 !!---------------------------------------------------------------------- 6 9 #if defined key_bdy 7 !!---------------------------------------------------------------------- -----------8 !! 'key_bdy' :Unstructured Open Boundary Conditions9 !!---------------------------------------------------------------------- -----------10 !!---------------------------------------------------------------------- 11 !! 'key_bdy' Unstructured Open Boundary Conditions 12 !!---------------------------------------------------------------------- 10 13 !! bdy_tra : Relaxation of tracers on unstructured open boundaries 11 !!--------------------------------------------------------------------------------- 12 !! * Modules used 14 !!---------------------------------------------------------------------- 13 15 USE oce ! ocean dynamics and tracers variables 14 16 USE dom_oce ! ocean space and time domain variables 15 17 USE bdy_oce ! ocean open boundary conditions 16 18 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 19 USE in_out_manager ! I/O manager 17 20 18 21 IMPLICIT NONE 19 22 PRIVATE 20 23 21 !! * Accessibility22 24 PUBLIC bdy_tra ! routine called in tranxt.F90 23 25 24 !! * Substitutions25 26 !! ---------------------------------------------------------------------------------27 !! OPA 9.0 , LODYC-IPSL (2003)28 !!---------------------------------------------------------------------- -----------26 !!---------------------------------------------------------------------- 27 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 28 !! $Id: $ 29 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 30 !!---------------------------------------------------------------------- 29 31 30 32 CONTAINS 31 33 32 34 SUBROUTINE bdy_tra( kt ) 33 !!---------------------------------------------------------------------- --------35 !!---------------------------------------------------------------------- 34 36 !! *** SUBROUTINE bdy_tra *** 35 37 !! … … 37 39 !! case of unstructured open boundaries. 38 40 !! 39 !! Reference : Engedahl H., 1995: Use of the flow relaxation scheme in 40 !! a three-dimensional baroclinic ocean model with realistic 41 !! topography. Tellus, 365-382. 42 !! History : 43 !! 9.0 ! 05-01 (J. Chanut, A. Sellar) Original 44 !!------------------------------------------------------------------------------ 45 !! * Arguments 41 !! Reference : Engedahl H., 1995, Tellus, 365-382. 42 !!---------------------------------------------------------------------- 46 43 INTEGER, INTENT( in ) :: kt 44 !! 45 REAL(wp) :: zwgt ! boundary weight 46 INTEGER :: ib, ik, igrd ! dummy loop indices 47 INTEGER :: ii, ij ! 2D addresses 48 !!---------------------------------------------------------------------- 49 ! 50 IF(ln_bdy_tra_frs) THEN ! If this is false, then this routine does nothing. 47 51 48 !! * Local declarations 49 REAL(wp) :: zwgt ! boundary weight 50 INTEGER :: jb, jk, jgrd ! dummy loop indices 51 INTEGER :: ii, ij ! 2D addresses 52 !!------------------------------------------------------------------------------ 53 54 jgrd=1 !: Everything is at T-points here 55 56 DO jb = 1, nblen(jgrd) 57 DO jk = 1, jpkm1 58 ii = nbi(jb,jgrd) 59 ij = nbj(jb,jgrd) 60 zwgt = nbw(jb,jgrd) 61 62 ! Temperature relaxation at the boundary 63 ta(ii,ij,jk) = ( ta(ii,ij,jk)*(1.-zwgt) + tbdy(jb,jk)*zwgt ) & 64 * tmask(ii,ij,jk) 65 66 ! Salinity relaxation at the boundary 67 sa(ii,ij,jk) = ( sa(ii,ij,jk)*(1.-zwgt) + sbdy(jb,jk)*zwgt ) & 68 * tmask(ii,ij,jk) 69 52 IF( kt == nit000 ) THEN 53 IF(lwp) WRITE(numout,*) 54 IF(lwp) WRITE(numout,*) 'bdy_tra : Flow Relaxation Scheme for tracers' 55 IF(lwp) WRITE(numout,*) '~~~~~~~' 56 ENDIF 57 ! 58 igrd = 1 ! Everything is at T-points here 59 DO ib = 1, nblen(igrd) 60 DO ik = 1, jpkm1 61 ii = nbi(ib,igrd) 62 ij = nbj(ib,igrd) 63 zwgt = nbw(ib,igrd) 64 ta(ii,ij,ik) = ( ta(ii,ij,ik) * (1.-zwgt) + tbdy(ib,ik) * zwgt ) * tmask(ii,ij,ik) 65 sa(ii,ij,ik) = ( sa(ii,ij,ik) * (1.-zwgt) + sbdy(ib,ik) * zwgt ) * tmask(ii,ij,ik) 70 66 END DO 71 67 END DO 72 73 CALL lbc_lnk( ta, 'T', 1. ) ! Boundary points should be updated 74 CALL lbc_lnk( sa, 'T', 1. ) ! 68 ! 69 CALL lbc_lnk( ta, 'T', 1. ) ! Boundary points should be updated 70 CALL lbc_lnk( sa, 'T', 1. ) ! 71 ! 72 ENDIF ! ln_bdy_tra_frs 75 73 76 74 END SUBROUTINE bdy_tra 75 77 76 #else 78 !!---------------------------------------------------------------------- -----------79 !! D efault option Empty module80 !!---------------------------------------------------------------------- -----------77 !!---------------------------------------------------------------------- 78 !! Dummy module NO Unstruct Open Boundary Conditions 79 !!---------------------------------------------------------------------- 81 80 CONTAINS 82 SUBROUTINE bdy_tra ! Empty routine 81 SUBROUTINE bdy_tra(kt) ! Empty routine 82 WRITE(*,*) 'bdy_tra: You should not have seen this print! error?', kt 83 83 END SUBROUTINE bdy_tra 84 84 #endif 85 85 86 !!====================================================================== ===========86 !!====================================================================== 87 87 END MODULE bdytra -
trunk/NEMO/OPA_SRC/BDY/bdyvol.F90
r911 r1125 1 1 MODULE bdyvol 2 !!====================================================================== ===========2 !!====================================================================== 3 3 !! *** MODULE bdyvol *** 4 4 !! Ocean dynamic : Volume constraint when unstructured boundary 5 5 !! and Free surface are used 6 !!================================================================================= 6 !!====================================================================== 7 !! History : 1.0 ! 2005-01 (J. Chanut, A. Sellar) Original code 8 !! - ! 2006-01 (J. Chanut) Bug correction 9 !! 3.0 ! 2008-04 (NEMO team) add in the reference version 10 !!---------------------------------------------------------------------- 7 11 #if defined key_bdy && defined key_dynspg_flt 8 !!--------------------------------------------------------------------------------- 9 !! 'key_bdy' and unstructured open boundary conditions 10 !! 'key_dynspg_flt' constant volume free surface 11 !!--------------------------------------------------------------------------------- 12 !! * Modules used 12 !!---------------------------------------------------------------------- 13 !! 'key_bdy' and unstructured open boundary conditions 14 !! 'key_dynspg_flt' filtered free surface 15 !!---------------------------------------------------------------------- 13 16 USE oce ! ocean dynamics and tracers 14 17 USE dom_oce ! ocean space and time domain … … 22 25 PRIVATE 23 26 24 !! * Accessibility25 27 PUBLIC bdy_vol ! routine called by dynspg_flt.h90 26 28 27 29 !! * Substitutions 28 30 # include "domzgr_substitute.h90" 29 !!--------------------------------------------------------------------------------- 30 !! OPA 9.0 , LODYC-IPSL (2003) 31 !!--------------------------------------------------------------------------------- 31 !!---------------------------------------------------------------------- 32 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 33 !! $Id: $ 34 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 35 !!---------------------------------------------------------------------- 32 36 33 37 CONTAINS 34 38 35 SUBROUTINE bdy_vol 36 !!---------------------------------------------------------------------- --------39 SUBROUTINE bdy_vol( kt ) 40 !!---------------------------------------------------------------------- 37 41 !! *** ROUTINE bdyvol *** 38 42 !! 39 !! ** Purpose : 40 !! This routine is called in dynspg_flt to control 43 !! ** Purpose : This routine is called in dynspg_flt to control 41 44 !! the volume of the system. A correction velocity is calculated 42 45 !! to correct the total transport through the unstructured OBC. … … 44 47 !! linear free surface coded in OPA 8.2 45 48 !! 46 !! ** Method : 47 !! The correction velocity (zubtpecor here) is defined calculating 49 !! ** Method : The correction velocity (zubtpecor here) is defined calculating 48 50 !! the total transport through all open boundaries (trans_bdy) minus 49 !! the cumulate E-P flux (z Cflxemp) divided by the total lateral51 !! the cumulate E-P flux (z_cflxemp) divided by the total lateral 50 52 !! surface (bdysurftot) of the unstructured boundary. 51 !! 52 !! zubtpecor = [trans_bdy - zCflxemp ]*(1./bdysurftot) 53 !! 54 !! with zCflxemp => sum of (Evaporation minus Precipitation) 53 !! zubtpecor = [trans_bdy - z_cflxemp ]*(1./bdysurftot) 54 !! with z_cflxemp => sum of (Evaporation minus Precipitation) 55 55 !! over all the domain in m3/s at each time step. 56 !! 57 !! zCflxemp < 0 when precipitation dominate 58 !! zCflxemp > 0 when evaporation dominate 56 !! z_cflxemp < 0 when precipitation dominate 57 !! z_cflxemp > 0 when evaporation dominate 59 58 !! 60 59 !! There are 2 options (user's desiderata): 61 !!62 60 !! 1/ The volume changes according to E-P, this is the default 63 61 !! option. In this case the cumulate E-P flux are setting to 64 !! zero (z Cflxemp=0) to calculate the correction velocity. So62 !! zero (z_cflxemp=0) to calculate the correction velocity. So 65 63 !! it will only balance the flux through open boundaries. 66 64 !! (set volbdy to 0 in tne namelist for this option) 67 !!68 65 !! 2/ The volume is constant even with E-P flux. In this case 69 66 !! the correction velocity must balance both the flux … … 71 68 !! surface. 72 69 !! (set volbdy to 1 in tne namelist for this option) 70 !!---------------------------------------------------------------------- 71 INTEGER, INTENT( in ) :: kt ! ocean time-step index 73 72 !! 74 !! History : 75 !! 8.5 ! 05-01 (J. Chanut, A. Sellar) Original code 76 !! ! 06-01 (J. Chanut) Bug correction 77 !!---------------------------------------------------------------------------- 78 !! * Arguments 79 INTEGER, INTENT( in ) :: kt ! ocean time-step index 80 81 !! * Local declarations 82 INTEGER :: ji,jj,jb, jgrd, jk 83 REAL(wp) :: zubtpecor 84 REAL(wp) :: zCflxemp 85 REAL(wp) :: ztranst 73 INTEGER :: ji, jj, jk, jb, jgrd 74 INTEGER :: ii, ij 75 REAL(wp) :: zubtpecor, z_cflxemp, ztranst 86 76 !!----------------------------------------------------------------------------- 87 77 88 78 IF( kt == nit000 ) THEN 89 IF(lwp) WRITE(numout,*) ' '79 IF(lwp) WRITE(numout,*) 90 80 IF(lwp) WRITE(numout,*)'bdy_vol : Correction of velocities along unstructured OBC' 91 81 IF(lwp) WRITE(numout,*)'~~~~~~~' 92 IF(lwp) WRITE(numout,*)' '93 82 END IF 94 83 95 ! 1. Calculate the cumulate surface Flux zCflxemp (m3/s) over all the domain. 96 ! --------------------------------------------------------------------------- 97 98 zCflxemp = 0.e0 99 84 ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain 85 ! ----------------------------------------------------------------------- 86 z_cflxemp = 0.e0 100 87 DO jj = 1, jpj 101 88 DO ji = 1, jpi 102 zCflxemp = zCflxemp + ( (emp(ji,jj)*bdytmask(ji,jj)*tmask_i(ji,jj) )/rauw) & 103 *e1v(ji,jj)*e2u(ji,jj) 89 z_cflxemp = z_cflxemp + emp(ji,jj) * bdytmask(ji,jj) * tmask_i(ji,jj) / rauw * e1v(ji,jj) * e2u(ji,jj) 104 90 END DO 105 91 END DO 106 IF( lk_mpp ) CALL mpp_sum( z Cflxemp ) ! sum over the global domain92 IF( lk_mpp ) CALL mpp_sum( z_cflxemp ) ! sum over the global domain 107 93 108 ! 2. Barotropic velocity through the unstructured open boundary 109 ! ------------------------------------------------------------- 110 94 ! Barotropic velocity through the unstructured open boundary 95 ! ---------------------------------------------------------- 111 96 zubtpecor = 0.e0 112 113 jgrd = 2 ! cumulate u component contribution first 97 jgrd = 2 ! cumulate u component contribution first 114 98 DO jb = 1, nblenrim(jgrd) 115 DO jk = 1, jpkm1116 zubtpecor = zubtpecor + flagu(jb) * ua (nbi(jb,jgrd), nbj(jb,jgrd), jk) &117 * e2u(nbi(jb,jgrd), nbj(jb,jgrd)) &118 * fse3u(nbi(jb,jgrd), nbj(jb,jgrd),jk)119 END DO99 DO jk = 1, jpkm1 100 ii = nbi(jb,jgrd) 101 ij = nbj(jb,jgrd) 102 zubtpecor = zubtpecor + flagu(jb) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 103 END DO 120 104 END DO 121 122 jgrd = 3 ! then add v component contribution 123 DO jb = 1, nblenrim(jgrd) 124 DO jk = 1, jpkm1 125 zubtpecor = zubtpecor + flagv(jb) * va (nbi(jb,jgrd), nbj(jb,jgrd), jk) & 126 * e1v(nbi(jb,jgrd), nbj(jb,jgrd)) & 127 * fse3v(nbi(jb,jgrd), nbj(jb,jgrd), jk) 128 END DO 105 jgrd = 3 ! then add v component contribution 106 DO jb = 1, nblenrim(jgrd) 107 DO jk = 1, jpkm1 108 ii = nbi(jb,jgrd) 109 ij = nbj(jb,jgrd) 110 zubtpecor = zubtpecor + flagv(jb) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 111 END DO 129 112 END DO 130 131 113 IF( lk_mpp ) CALL mpp_sum( zubtpecor ) ! sum over the global domain 132 114 133 134 ! 3. The normal velocity correction 135 ! --------------------------------- 136 137 IF (volbdy==1) THEN 138 zubtpecor = (zubtpecor - zCflxemp)*(1./bdysurftot) 139 ELSE 140 zubtpecor = zubtpecor*(1./bdysurftot) 115 ! The normal velocity correction 116 ! ------------------------------ 117 IF (volbdy==1) THEN ; zubtpecor = ( zubtpecor - z_cflxemp) / bdysurftot 118 ELSE ; zubtpecor = zubtpecor / bdysurftot 141 119 END IF 142 120 121 ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation 122 ! ------------------------------------------------------------- 123 ztranst = 0.e0 124 jgrd = 2 ! correct u component 125 DO jb = 1, nblenrim(jgrd) 126 DO jk = 1, jpkm1 127 ii = nbi(jb,jgrd) 128 ij = nbj(jb,jgrd) 129 ua(ii,ij,jk) = ua(ii,ij,jk) - flagu(jb) * zubtpecor * umask(ii,ij,jk) 130 ztranst = ztranst + flagu(jb) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 131 END DO 132 END DO 133 jgrd = 3 ! correct v component 134 DO jb = 1, nblenrim(jgrd) 135 DO jk = 1, jpkm1 136 ii = nbi(jb,jgrd) 137 ij = nbj(jb,jgrd) 138 va(ii,ij,jk) = va(ii,ij,jk) -flagv(jb) * zubtpecor * vmask(ii,ij,jk) 139 ztranst = ztranst + flagv(jb) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 140 END DO 141 END DO 142 IF( lk_mpp ) CALL mpp_sum( ztranst ) ! sum over the global domain 143 144 ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected 145 ! ------------------------------------------------------ 146 143 147 IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 144 IF(lwp) WRITE(numout,*) ' '148 IF(lwp) WRITE(numout,*) 145 149 IF(lwp) WRITE(numout,*)'bdy_vol : time step :', kt 146 150 IF(lwp) WRITE(numout,*)'~~~~~~~ ' 147 IF(lwp) WRITE(numout,*)' ' 148 IF(lwp) WRITE(numout,*)' cumulate flux EMP :', zCflxemp,' (m3/s)' 149 IF(lwp) WRITE(numout,*)' total lateral surface of OBC :',bdysurftot,'(m2)' 150 IF(lwp) WRITE(numout,*)' correction velocity zubtpecor :',zubtpecor,'(m/s)' 151 IF(lwp) WRITE(numout,*)' ' 151 IF(lwp) WRITE(numout,*)' cumulate flux EMP =', z_cflxemp , ' (m3/s)' 152 IF(lwp) WRITE(numout,*)' total lateral surface of OBC =', bdysurftot, '(m2)' 153 IF(lwp) WRITE(numout,*)' correction velocity zubtpecor =', zubtpecor , '(m/s)' 154 IF(lwp) WRITE(numout,*)' cumulated transport ztranst =', ztranst , '(m3/s)' 152 155 END IF 153 154 ! 4. Correction of the total velocity on the unstructured 155 ! boundary to respect the mass flux conservation 156 ! ------------------------------------------------------- 157 158 ztranst = 0.e0 159 160 jgrd = 2 ! correct u component 161 DO jb = 1, nblenrim(jgrd) 162 DO jk = 1, jpkm1 163 ua(nbi(jb, jgrd), nbj(jb, jgrd), jk) = ua(nbi(jb, jgrd), nbj(jb, jgrd), jk) & 164 -flagu(jb) * zubtpecor * umask(nbi(jb, jgrd), nbj(jb, jgrd), jk) 165 ztranst = ztranst + flagu(jb) * ua (nbi(jb,jgrd), nbj(jb,jgrd), jk) & 166 * e2u(nbi(jb,jgrd), nbj(jb,jgrd)) & 167 * fse3u(nbi(jb,jgrd), nbj(jb,jgrd), jk) 168 END DO 169 END DO 170 171 jgrd = 3 ! correct v component 172 DO jb = 1, nblenrim(jgrd) 173 DO jk = 1, jpkm1 174 va(nbi(jb, jgrd), nbj(jb, jgrd), jk) = va(nbi(jb, jgrd), nbj(jb, jgrd), jk) & 175 -flagv(jb) * zubtpecor * vmask(nbi(jb, jgrd), nbj(jb, jgrd), jk) 176 ztranst = ztranst + flagv(jb) * va (nbi(jb,jgrd), nbj(jb,jgrd), jk) & 177 * e1v(nbi(jb,jgrd), nbj(jb,jgrd)) & 178 * fse3v(nbi(jb,jgrd), nbj(jb,jgrd), jk) 179 END DO 180 END DO 181 182 IF( lk_mpp ) CALL mpp_sum( ztranst ) ! sum over the global domain 183 184 ! 5. Check the cumulate transport through unstructured OBC 185 ! once barotropic velocities corrected 186 ! -------------------------------------------------------- 187 188 IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 189 IF(lwp) WRITE(numout,*)' ' 190 IF(lwp) WRITE(numout,*)' Cumulate transport ztranst =', ztranst,'(m3/s)' 191 IF(lwp) WRITE(numout,*)' ' 192 END IF 193 156 ! 194 157 END SUBROUTINE bdy_vol 195 158 196 159 #else 197 !!---------------------------------------------------------------------- -----------198 !! Default option : Empty module199 !!---------------------------------------------------------------------- -----------160 !!---------------------------------------------------------------------- 161 !! Dummy module NO Unstruct Open Boundary Conditions 162 !!---------------------------------------------------------------------- 200 163 CONTAINS 201 SUBROUTINE bdy_vol ! Empty routine 164 SUBROUTINE bdy_vol( kt ) ! Empty routine 165 WRITE(*,*) 'bdy_vol: You should not have seen this print! error?', kt 202 166 END SUBROUTINE bdy_vol 203 167 #endif 204 168 205 !!====================================================================== ===========169 !!====================================================================== 206 170 END MODULE bdyvol -
trunk/NEMO/OPA_SRC/DOM/domvvl.F90
r1058 r1125 302 302 #endif 303 303 304 #if defined key_bdy || defined key_bdy_tides 305 DO jj = 1, jpj 306 DO ji = 1, jpi 307 zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj) 308 END DO 309 END DO 304 #if defined key_bdy 305 zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:) 310 306 #endif 311 307 -
trunk/NEMO/OPA_SRC/DYN/divcur.F90
r1058 r1125 136 136 #endif 137 137 #endif 138 #if defined key_bdy || defined key_bdy_tides138 #if defined key_bdy 139 139 ! unstructured open boundaries (div must be zero behind the open boundary) 140 140 DO jj = 1, jpj … … 354 354 #endif 355 355 #endif 356 #if defined key_bdy || defined key_bdy_tides356 #if defined key_bdy 357 357 ! unstructured open boundaries (div must be zero behind the open boundary) 358 358 DO jj = 1, jpj -
trunk/NEMO/OPA_SRC/DYN/dynnxt.F90
r911 r1125 238 238 DO jk = 1, jpkm1 ! Horizontal slab 239 239 ! ! =============== 240 # elif defined key_bdy || key_bdy_tides240 # elif defined key_bdy 241 241 ! ! =============== 242 242 END DO ! End of slab -
trunk/NEMO/OPA_SRC/DYN/wzvmod.F90
r1040 r1125 61 61 REAL(wp) :: z2dt, zraur ! temporary scalar 62 62 REAL(wp), DIMENSION (jpi,jpj) :: zssha, zun, zvn, zhdiv 63 #if defined key_bdy || defined key_bdy_tides63 #if defined key_bdy 64 64 INTEGER :: jgrd, jb ! temporary scalars 65 65 #endif … … 113 113 #endif 114 114 115 #if defined key_bdy || defined key_bdy_tides115 #if defined key_bdy 116 116 jgrd=1 !: tracer grid. 117 117 DO jb = 1, nblenrim(jgrd) -
trunk/NEMO/OPA_SRC/opa.F90
r1037 r1125 53 53 USE bdy_par ! unstructured open boundary cond. parameters 54 54 USE bdyini ! unstructured open boundary cond. initialization (bdy_init routine) 55 USE bdytides ! tides at open boundaries initialization (lk_bdy_tides)56 55 USE istate ! initial state setting (istate_init routine) 57 56 USE eosbn2 ! equation of state (eos bn2 routine) … … 272 271 IF( lk_obc ) CALL obc_init ! Open boundaries 273 272 274 IF( lk_bdy .OR. lk_bdy_tides ) & 275 & CALL bdy_init ! Unstructured open boundaries 273 IF( lk_bdy ) CALL bdy_init ! Unstructured open boundaries 276 274 277 275 CALL istate_init ! ocean initial state (Dynamics and tracers)
Note: See TracChangeset
for help on using the changeset viewer.
