Changeset 2814 for branches/2011/UKMO_MERCATOR_obc_bdy_merge
- Timestamp:
- 2011-07-27T14:41:28+02:00 (13 years ago)
- Location:
- branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r2800 r2814 26 26 USE zdfbfr ! bottom friction 27 27 USE dynvor ! vorticity term 28 USE obc_par ! for lk_obc 28 29 USE obc_oce ! Lateral open boundary condition 29 30 USE obcdta ! open boundary condition data … … 528 529 ! ! ==================== ! 529 530 530 #if defined key_obc531 IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) !!gm totally useless ?????532 IF( lp_obc_west ) sshfow_b(:,:) = zcoef * sshfow_b(:,:)533 IF( lp_obc_north ) sshfon_b(:,:) = zcoef * sshfon_b(:,:)534 IF( lp_obc_south ) sshfos_b(:,:) = zcoef * sshfos_b(:,:)535 #endif536 537 531 ! ----------------------------------------------------------------------------- 538 532 ! Phase 3. update the general trend with the barotropic trend -
branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/FLO/floblk.F90
r2715 r2814 345 345 ! more time. 346 346 # if defined key_obc 347 DO jfl = 1, jpnfl 348 IF( lp_obc_east ) THEN 349 IF( jped <= zgjfl(jfl) .AND. zgjfl(jfl) <= jpef .AND. nieob-1 <= zgifl(jfl) ) THEN 350 zgifl (jfl) = INT(zgifl(jfl)) + 0.5 351 zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 352 zagefl(jfl) = rdt 353 END IF 354 END IF 355 IF( lp_obc_west ) THEN 356 IF( jpwd <= zgjfl(jfl) .AND. zgjfl(jfl) <= jpwf .AND. niwob >= zgifl(jfl) ) THEN 357 zgifl (jfl) = INT(zgifl(jfl)) + 0.5 358 zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 359 zagefl(jfl) = rdt 360 END IF 361 END IF 362 IF( lp_obc_north ) THEN 363 IF( jpnd <= zgifl(jfl) .AND. zgifl(jfl) <= jpnf .AND. njnob-1 >= zgjfl(jfl) ) THEN 364 zgifl (jfl) = INT(zgifl(jfl)) + 0.5 365 zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 366 zagefl(jfl) = rdt 367 END IF 368 END IF 369 IF( lp_obc_south ) THEN 370 IF( jpsd <= zgifl(jfl) .AND. zgifl(jfl) <= jpsf .AND. njsob >= zgjfl(jfl) ) THEN 371 zgifl (jfl) = INT(zgifl(jfl)) + 0.5 372 zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 373 zagefl(jfl) = rdt 374 END IF 375 END IF 376 END DO 347 !!!!!!!! NEED TO SORT THIS OUT !!!!!!!! 348 !!$ DO jfl = 1, jpnfl 349 !!$ IF( lp_obc_east ) THEN 350 !!$ IF( jped <= zgjfl(jfl) .AND. zgjfl(jfl) <= jpef .AND. nieob-1 <= zgifl(jfl) ) THEN 351 !!$ zgifl (jfl) = INT(zgifl(jfl)) + 0.5 352 !!$ zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 353 !!$ zagefl(jfl) = rdt 354 !!$ END IF 355 !!$ END IF 356 !!$ IF( lp_obc_west ) THEN 357 !!$ IF( jpwd <= zgjfl(jfl) .AND. zgjfl(jfl) <= jpwf .AND. niwob >= zgifl(jfl) ) THEN 358 !!$ zgifl (jfl) = INT(zgifl(jfl)) + 0.5 359 !!$ zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 360 !!$ zagefl(jfl) = rdt 361 !!$ END IF 362 !!$ END IF 363 !!$ IF( lp_obc_north ) THEN 364 !!$ IF( jpnd <= zgifl(jfl) .AND. zgifl(jfl) <= jpnf .AND. njnob-1 >= zgjfl(jfl) ) THEN 365 !!$ zgifl (jfl) = INT(zgifl(jfl)) + 0.5 366 !!$ zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 367 !!$ zagefl(jfl) = rdt 368 !!$ END IF 369 !!$ END IF 370 !!$ IF( lp_obc_south ) THEN 371 !!$ IF( jpsd <= zgifl(jfl) .AND. zgifl(jfl) <= jpsf .AND. njsob >= zgjfl(jfl) ) THEN 372 !!$ zgifl (jfl) = INT(zgifl(jfl)) + 0.5 373 !!$ zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 374 !!$ zagefl(jfl) = rdt 375 !!$ END IF 376 !!$ END IF 377 !!$ END DO 377 378 #endif 378 379 -
branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obc_oce.F90
r2800 r2814 56 56 ! !: =F read obc coordinates from namelist 57 57 LOGICAL :: ln_mask_file !: =T read obcmask from file 58 LOGICAL, DIMENSION(jp_obc) :: ln_tides !: =T apply tidal harmonic forcing along open boundaries59 58 LOGICAL :: ln_vol !: =T volume correction 60 59 LOGICAL, DIMENSION(jp_obc) :: ln_clim !: =T obc data files contain climatological data (time-cyclic) … … 64 63 INTEGER, DIMENSION(jp_obc) :: nn_dtactl !: = 0 use the initial state as obc dta ; 65 64 !: = 1 read it in a NetCDF file 65 INTEGER, DIMENSION(jp_obc) :: nn_tides !: = 0 no tidal harmonic forcing 66 !: = 1 apply ONLY tidal harmonic forcing for barotropic solution 67 !: = 2 ADD tidal harmonic forcing to other barotropic boundary data 66 68 INTEGER :: nn_volctl !: = 0 the total volume will have the variability of the surface Flux E-P 67 69 ! ! = 1 the volume will be constant during all the integration. … … 98 100 !!---------------------------------------------------------------------- 99 101 100 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) 101 TYPE(OBC_INDEX), DIMENSION(jp_obc), TARGET :: idx_obc !: obc indices (local process)102 TYPE(OBC_DATA) , DIMENSION(jp_obc) :: dta_obc !: obc external data (local process)102 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global !: workspace for reading in global data arrays 103 TYPE(OBC_INDEX), DIMENSION(jp_obc), TARGET :: idx_obc !: obc indices (local process) 104 TYPE(OBC_DATA) , DIMENSION(jp_obc) :: dta_obc !: obc external data (local process) 103 105 104 106 !!---------------------------------------------------------------------- -
branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdta.F90
r2800 r2814 65 65 !!---------------------------------------------------------------------- 66 66 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 67 USE wrk_nemo, ONLY: wrk_2d_ 9, wrk_2d_10! 2D workspace67 USE wrk_nemo, ONLY: wrk_2d_22, wrk_2d_23 ! 2D workspace 68 68 !! 69 69 INTEGER, INTENT( in ) :: kt ! ocean time-step index … … 76 76 !!--------------------------------------------------------------------------- 77 77 78 IF(wrk_in_use(2, 9,10) ) THEN78 IF(wrk_in_use(2, 22,23) ) THEN 79 79 CALL ctl_stop('obc_dta: ERROR: requested workspace arrays are unavailable.') ; RETURN 80 80 END IF … … 87 87 ! Calculate depth-mean currents 88 88 !----------------------------- 89 pu2d => wrk_2d_ 990 pu2d => wrk_2d_ 1089 pu2d => wrk_2d_22 90 pu2d => wrk_2d_23 91 91 92 92 pu2d(:,:) = 0.e0 … … 206 206 ! Update barotropic boundary conditions only 207 207 ! jit is optional argument for fld_read 208 IF( nn_dyn2d(ib_obc) .gt. 0 ) THEN208 IF( nn_dyn2d(ib_obc) .gt. 0 .and. nn_tides(ib_obc) .ne. 1 ) THEN 209 209 jend = jstart + 2 210 210 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), jit=jit ) 211 211 ENDIF 212 IF( nn_tides(ib_obc) .gt. 0 ) THEN 213 CALL tide_update( kt=kt, jit=jit, idx=idx_obc(ib_obc), dta=dta_obc(ib_obc), td=tides(ib_obc) ) 214 ENDIF 212 215 ELSE 213 jend = jstart + nb_obc_fld(ib_obc) - 1 214 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend ), map=nbmap_ptr(jstart:jend), timeshift=1 ) 216 IF( nb_obc_fld(ib_obc) .gt. 0 ) THEN 217 jend = jstart + nb_obc_fld(ib_obc) - 1 218 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), timeshift=1 ) 219 ENDIF 220 IF( nn_tides(ib_obc) .gt. 0 ) THEN 221 !!! THINK ABOUT kt, jit VALUES !!! 222 CALL tide_update( kt=kt, jit=0, idx=idx_obc(ib_obc), dta=dta_obc(ib_obc), td=tides(ib_obc) ) 223 ENDIF 215 224 ENDIF 216 225 jstart = jend+1 … … 257 266 END DO ! ib_obc 258 267 259 IF(wrk_not_released(2, 9,10) ) CALL ctl_stop('obc_dta: ERROR: failed to release workspace arrays.')268 IF(wrk_not_released(2, 22,23) ) CALL ctl_stop('obc_dta: ERROR: failed to release workspace arrays.') 260 269 261 270 END SUBROUTINE obc_dta … … 383 392 IF( nn_dyn2d(ib_obc) .gt. 0 ) THEN 384 393 385 IF( nn_dyn2d(ib_obc) .ne. jp_frs ) THEN394 IF( nn_dyn2d(ib_obc) .ne. jp_frs .and. nn_tides(ib_obc) .ne. 1) THEN 386 395 jfld = jfld + 1 387 396 blf_i(jfld) = bn_ssh … … 392 401 ENDIF 393 402 394 IF( .not. ln_full_vel_array(ib_obc) ) THEN403 IF( .not. ln_full_vel_array(ib_obc) .and. nn_tides(ib_obc) .ne. 1 ) THEN 395 404 396 405 jfld = jfld + 1 … … 604 613 605 614 ! nn_dtactl = 1 606 ! Set boundary data arrays to point to relevant"fnow" arrays607 !--------------------------------------------------- --------615 ! Set boundary data arrays to point to "fnow" arrays 616 !--------------------------------------------------- 608 617 IF (nn_dyn2d(ib_obc) .gt. 0) THEN 609 IF( nn_dyn2d(ib_obc) .ne. jp_frs ) THEN618 IF( nn_dyn2d(ib_obc) .ne. jp_frs .and. nn_tides(ib_obc) .ne. 1 ) THEN 610 619 jfld = jfld + 1 611 620 dta_obc(ib_obc)%ssh => bf(jfld)%fnow(:,1,1) 612 621 ENDIF 613 IF( ln_full_vel_array(ib_obc) ) THEN622 IF( ln_full_vel_array(ib_obc) .or. nn_tides(ib_obc) .eq. 1 ) THEN 614 623 ! In this case we need space but we aren't reading it 615 624 ! directly from the external file. 616 625 IF( nn_dyn2d(ib_obc) .eq. jp_frs ) THEN 617 ilen0(2) = nblen(2) 618 ilen0(3) = nblen(3) 626 ilen0(:) = nblen(:) 619 627 ELSE 620 ilen0(2) = nblenrim(2) 621 ilen0(3) = nblenrim(3) 628 ilen0(:) = nblenrim(:) 622 629 ENDIF 630 IF( nn_tides(ib_obc) .eq. 1 ) ALLOCATE( dta_obc(ib_obc)%ssh(ilen0(1)) ) 623 631 ALLOCATE( dta_obc(ib_obc)%u2d(ilen0(2)) ) 624 632 ALLOCATE( dta_obc(ib_obc)%v2d(ilen0(3)) ) … … 631 639 ENDIF 632 640 IF (nn_dyn3d(ib_obc) .gt. 0 .or. & 633 & ln_full_vel_array(ib_obc) .and. nn_dyn2d(ib_obc) .gt. 0) THEN641 & ( ln_full_vel_array(ib_obc) .and. nn_dyn2d(ib_obc) .gt. 0 ) ) THEN 634 642 jfld = jfld + 1 635 643 dta_obc(ib_obc)%u3d => bf(jfld)%fnow(:,1,:) … … 666 674 CONTAINS 667 675 SUBROUTINE obc_dta( kt, jit ) ! Empty routine 676 INTEGER, INTENT( in ) :: kt 677 INTEGER, INTENT( in ), OPTIONAL :: jit 668 678 WRITE(*,*) 'obc_dta: You should not have seen this print! error?', kt 669 679 END SUBROUTINE obc_dta -
branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn2d.F90
r2800 r2814 24 24 USE phycst ! physical constants 25 25 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 26 USE obctides ! for tidal harmonic forcing at boundary27 26 USE in_out_manager ! 28 27 … … 30 29 PRIVATE 31 30 32 PUBLIC obc_dyn2d ! routine called in dynspg_ts (time splitting case ONLY)31 PUBLIC obc_dyn2d ! routine called in dynspg_ts and dyn_nxt and dynspg_flt 33 32 34 33 !!---------------------------------------------------------------------- … … 77 76 !! topography. Tellus, 365-382. 78 77 !!---------------------------------------------------------------------- 79 INTEGER, INTENT(in) :: kt78 INTEGER, INTENT(in) :: kt 80 79 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 81 80 TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data … … 114 113 !! 115 114 !! - Apply Flather boundary conditions on normal barotropic velocities 116 !! (ln_dyn_fla=.true. or ln_tides=.true.)117 115 !! 118 116 !! ** WARNINGS about FLATHER implementation: … … 143 141 ! ---------------------------------! 144 142 145 !!!!!!!!!!!! SOME WORK TO DO ON THE TIDES HERE !!!!!!!!!!!!!146 147 143 !!! REPLACE spgu with nemo_wrk work space 148 144 … … 154 150 ij = idx%nbj(jb,igrd) 155 151 spgu(ii, ij) = dta%ssh(jb) 156 !!$ IF( ln_tides ) spgu(ii, ij) = spgu(ii, ij) + sshtide(jb)157 152 END DO 158 153 ! … … 167 162 ! 168 163 zcorr = - idx%flagu(jb) * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 169 !!$ zforc = dta%u2d(jb) + utide(jb)170 164 zforc = dta%u2d(jb) 171 165 pu2d(ii,ij) = zforc + zcorr * umask(ii,ij,1) … … 181 175 ! 182 176 zcorr = - idx%flagv(jb) * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 183 !!$ zforc = dta%v2d(jb) + vtide(jb)184 177 zforc = dta%v2d(jb) 185 178 pv2d(ii,ij) = zforc + zcorr * vmask(ii,ij,1) -
branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcini.F90
r2797 r2814 22 22 USE dom_oce ! ocean space and time domain 23 23 USE obc_oce ! unstructured open boundary conditions 24 USE obctides ! tides at open boundaries initialization (tide_init routine)25 24 USE in_out_manager ! I/O units 26 25 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 73 72 & nn_ice_lim2, & 74 73 #endif 75 & ln_tides, ln_vol, ln_clim, nn_dtactl, nn_volctl, &74 & nn_tides, ln_vol, ln_clim, nn_dtactl, nn_volctl, & 76 75 & nn_rimwidth, nn_dmp2d_in, nn_dmp2d_out, & 77 76 & nn_dmp3d_in, nn_dmp3d_out … … 105 104 nn_ice_lim2(:) = 0 106 105 #endif 107 ln_tides(:) = .false.108 106 ln_vol = .false. 109 107 ln_clim(:) = .false. 110 108 nn_dtactl(:) = -1 ! uninitialised flag 109 nn_tides(:) = 0 ! default to no tidal forcing 111 110 nn_volctl = -1 ! uninitialised flag 112 111 nn_rimwidth(:) = -1 ! uninitialised flag … … 182 181 IF(lwp) WRITE(numout,*) 183 182 184 IF( ln_tides(ib_obc) ) THEN 185 IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing at unstructured open boundaries' 183 SELECT CASE( nn_tides(ib_obc) ) 184 CASE(0) 185 IF(lwp) WRITE(numout,*) 'No tidal harmonic forcing' 186 186 IF(lwp) WRITE(numout,*) 187 ENDIF 188 189 !!$ ! Presumably will need to read in a separate namelist for each boundary that includes tides??? 190 !!$ IF( ln_tides ) CALL tide_init( ib_obc ) ! Read tides namelist 191 187 CASE(1) 188 IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing ONLY for barotropic solution' 189 IF(lwp) WRITE(numout,*) 190 CASE(2) 191 IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing ADDED to other barotropic boundary conditions' 192 IF(lwp) WRITE(numout,*) 193 CASE DEFAULT 194 CALL ctl_stop( 'obc_ini: ERROR: incorrect value for nn_tides ' ) 195 END SELECT 192 196 193 197 ENDDO … … 345 349 ENDDO 346 350 347 ! Compute rim weights 348 ! ------------------- 351 ! Compute rim weights for FRS scheme 352 ! ---------------------------------- 349 353 DO igrd = 1, jpbgrd 350 354 DO ib = 1, idx_obc(ib_obc)%nblen(igrd) … … 537 541 IF( lk_mpp ) CALL mpp_sum( obcsurftot ) ! sum over the global domain 538 542 END IF 539 540 ! Read in tidal constituents and adjust for model start time541 ! ----------------------------------------------------------542 !!$ IF( ln_tides ) CALL tide_data543 543 ! 544 544 ! Tidy up … … 550 550 #else 551 551 !!--------------------------------------------------------------------------------- 552 !! Dummy module NO unstructuredopen boundaries552 !! Dummy module NO open boundaries 553 553 !!--------------------------------------------------------------------------------- 554 554 CONTAINS -
branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obctides.F90
r2797 r2814 10 10 !!---------------------------------------------------------------------- 11 11 #if defined key_obc 12 !!$ !!---------------------------------------------------------------------- 13 !!$ !! 'key_obc' Unstructured Open Boundary Condition 14 !!$ !!---------------------------------------------------------------------- 15 !!$ !! PUBLIC 16 !!$ !! tide_init : read of namelist 17 !!$ !! tide_data : read in and initialisation of tidal constituents at boundary 18 !!$ !! tide_update : calculation of tidal forcing at each timestep 19 !!$ !! PRIVATE 20 !!$ !! uvset :\ 21 !!$ !! vday : | Routines to correct tidal harmonics forcing for 22 !!$ !! shpen : | start time of integration 23 !!$ !! ufset : | 24 !!$ !! vset :/ 25 !!$ !!---------------------------------------------------------------------- 26 !!$ USE oce ! ocean dynamics and tracers 27 !!$ USE dom_oce ! ocean space and time domain 28 !!$ USE iom 29 !!$ USE in_out_manager ! I/O units 30 !!$ USE phycst ! physical constants 31 !!$ USE lbclnk ! ocean lateral boundary conditions (or mpp link) 32 !!$ USE obc_par ! Unstructured boundary parameters 33 !!$ USE obc_oce ! ocean open boundary conditions 34 !!$ USE daymod ! calendar 35 !!$ 36 !!$ IMPLICIT NONE 37 !!$ PRIVATE 38 !!$ 39 !!$ PUBLIC tide_init ! routine called in obcini 40 !!$ PUBLIC tide_data ! routine called in obcini 41 !!$ PUBLIC tide_update ! routine called in obcdyn 42 !!$ 43 !!$ LOGICAL, PUBLIC :: ln_tide_date !: =T correct tide phases and amplitude for model start date 44 !!$ INTEGER, PUBLIC, PARAMETER :: jptides_max = 15 !: Max number of tidal contituents 45 !!$ INTEGER, PUBLIC :: 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 , PUBLIC, DIMENSION(jptides_max) :: nindx !: ??? 51 !!$ REAL(wp), PUBLIC, DIMENSION(jptides_max) :: tide_speed !: Phase speed of tidal constituent (deg/hr) 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 56 !!$ 57 !!$ !!---------------------------------------------------------------------- 58 !!$ !! NEMO/OPA 3.3 , NEMO Consortium (2010) 59 !!$ !! $Id: obctides.F90 2528 2010-12-27 17:33:53Z rblod $ 60 !!$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 61 !!$ !!---------------------------------------------------------------------- 62 !!$CONTAINS 63 !!$ 64 !!$ SUBROUTINE tide_init 65 !!$ !!---------------------------------------------------------------------- 66 !!$ !! *** SUBROUTINE tide_init *** 67 !!$ !! 68 !!$ !! ** Purpose : - Read in namelist for tides 69 !!$ !! 70 !!$ !!---------------------------------------------------------------------- 71 !!$ INTEGER :: itide ! dummy loop index 72 !!$ !! 73 !!$ NAMELIST/namobc_tide/ln_tide_date, filtide, tide_cpt, tide_speed 74 !!$ !!---------------------------------------------------------------------- 75 !!$ 76 !!$ IF(lwp) WRITE(numout,*) 77 !!$ IF(lwp) WRITE(numout,*) 'tide_init : initialization of tidal harmonic forcing at open boundaries' 78 !!$ IF(lwp) WRITE(numout,*) '~~~~~~~~~' 79 !!$ 80 !!$ ! Namelist namobc_tide : tidal harmonic forcing at open boundaries 81 !!$ ln_tide_date = .false. 82 !!$ filtide(:) = '' 83 !!$ tide_speed(:) = 0.0 84 !!$ tide_cpt(:) = '' 85 !!$ REWIND( numnam ) ! Read namelist parameters 86 !!$ READ ( numnam, namobc_tide ) 87 !!$ ! ! Count number of components specified 88 !!$ ntide = jptides_max 89 !!$ DO itide = 1, jptides_max 90 !!$ IF( tide_cpt(itide) == '' ) THEN 91 !!$ ntide = itide-1 92 !!$ exit 93 !!$ ENDIF 94 !!$ END DO 95 !!$ 96 !!$ ! ! find constituents in standard list 97 !!$ DO itide = 1, ntide 98 !!$ nindx(itide) = 0 99 !!$ IF( TRIM( tide_cpt(itide) ) == 'Q1' ) nindx(itide) = 1 100 !!$ IF( TRIM( tide_cpt(itide) ) == 'O1' ) nindx(itide) = 2 101 !!$ IF( TRIM( tide_cpt(itide) ) == 'P1' ) nindx(itide) = 3 102 !!$ IF( TRIM( tide_cpt(itide) ) == 'S1' ) nindx(itide) = 4 103 !!$ IF( TRIM( tide_cpt(itide) ) == 'K1' ) nindx(itide) = 5 104 !!$ IF( TRIM( tide_cpt(itide) ) == '2N2' ) nindx(itide) = 6 105 !!$ IF( TRIM( tide_cpt(itide) ) == 'MU2' ) nindx(itide) = 7 106 !!$ IF( TRIM( tide_cpt(itide) ) == 'N2' ) nindx(itide) = 8 107 !!$ IF( TRIM( tide_cpt(itide) ) == 'NU2' ) nindx(itide) = 9 108 !!$ IF( TRIM( tide_cpt(itide) ) == 'M2' ) nindx(itide) = 10 109 !!$ IF( TRIM( tide_cpt(itide) ) == 'L2' ) nindx(itide) = 11 110 !!$ IF( TRIM( tide_cpt(itide) ) == 'T2' ) nindx(itide) = 12 111 !!$ IF( TRIM( tide_cpt(itide) ) == 'S2' ) nindx(itide) = 13 112 !!$ IF( TRIM( tide_cpt(itide) ) == 'K2' ) nindx(itide) = 14 113 !!$ IF( TRIM( tide_cpt(itide) ) == 'M4' ) nindx(itide) = 15 114 !!$ IF( nindx(itide) == 0 .AND. lwp ) THEN 115 !!$ WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list' 116 !!$ CALL ctl_warn( ctmp1 ) 117 !!$ ENDIF 118 !!$ END DO 119 !!$ ! ! Parameter control and print 120 !!$ IF( ntide < 1 ) THEN 121 !!$ CALL ctl_stop( ' Did not find any tidal components in namelist namobc_tide' ) 122 !!$ ELSE 123 !!$ IF(lwp) WRITE(numout,*) ' Namelist namobc_tide : tidal harmonic forcing at open boundaries' 124 !!$ IF(lwp) WRITE(numout,*) ' tidal components specified ', ntide 125 !!$ IF(lwp) WRITE(numout,*) ' ', tide_cpt(1:ntide) 126 !!$ IF(lwp) WRITE(numout,*) ' associated phase speeds (deg/hr) : ' 127 !!$ IF(lwp) WRITE(numout,*) ' ', tide_speed(1:ntide) 128 !!$ ENDIF 129 !!$ 130 !!$ ! Initialisation of tidal harmonics arrays 131 !!$ sshtide(:) = 0.e0 132 !!$ utide (:) = 0.e0 133 !!$ vtide (:) = 0.e0 134 !!$ ! 135 !!$ END SUBROUTINE tide_init 136 !!$ 137 !!$ 138 !!$ SUBROUTINE tide_data 139 !!$ !!---------------------------------------------------------------------- 140 !!$ !! *** SUBROUTINE tide_data *** 141 !!$ !! 142 !!$ !! ** Purpose : - Read in tidal harmonics data and adjust for the start 143 !!$ !! time of the model run. 144 !!$ !! 145 !!$ !!---------------------------------------------------------------------- 146 !!$ INTEGER :: itide, igrd, ib ! dummy loop indices 147 !!$ CHARACTER(len=80) :: clfile ! full file name for tidal input file 148 !!$ INTEGER :: ipi, ipj, inum, idvar ! temporary integers (netcdf read) 149 !!$ INTEGER, DIMENSION(6) :: lendta=0 ! length of data in the file (note may be different from nblendta!) 150 !!$ REAL(wp) :: z_arg, z_atde, z_btde, z1t, z2t 151 !!$ REAL(wp), DIMENSION(jpbdta,1) :: zdta ! temporary array for data fields 152 !!$ REAL(wp), DIMENSION(jptides_max) :: z_vplu, z_ftc 153 !!$ !!------------------------------------------------------------------------------ 154 !!$ 155 !!$ ! Open files and read in tidal forcing data 156 !!$ ! ----------------------------------------- 157 !!$ 158 !!$ ipj = 1 159 !!$ 160 !!$ DO itide = 1, ntide 161 !!$ ! ! SSH fields 162 !!$ clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc' 163 !!$ IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 164 !!$ CALL iom_open( clfile, inum ) 165 !!$ igrd = 4 166 !!$ IF( nblendta(igrd) <= 0 ) THEN 167 !!$ idvar = iom_varid( inum,'z1' ) 168 !!$ IF(lwp) WRITE(numout,*) 'iom_file(1)%ndims(idvar) : ',iom_file%ndims(idvar) 169 !!$ nblendta(igrd) = iom_file(inum)%dimsz(1,idvar) 170 !!$ WRITE(numout,*) 'Dim size for z1 is ', nblendta(igrd) 171 !!$ ENDIF 172 !!$ ipi = nblendta(igrd) 173 !!$ CALL iom_get( inum, jpdom_unknown, 'z1', zdta(1:ipi,1:ipj) ) 174 !!$ DO ib = 1, nblenrim(igrd) 175 !!$ ssh1(ib,itide) = zdta(nbmap(ib,igrd),1) 176 !!$ END DO 177 !!$ CALL iom_get( inum, jpdom_unknown, 'z2', zdta(1:ipi,1:ipj) ) 178 !!$ DO ib = 1, nblenrim(igrd) 179 !!$ ssh2(ib,itide) = zdta(nbmap(ib,igrd),1) 180 !!$ END DO 181 !!$ CALL iom_close( inum ) 182 !!$ ! 183 !!$ ! ! U fields 184 !!$ clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc' 185 !!$ IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 186 !!$ CALL iom_open( clfile, inum ) 187 !!$ igrd = 5 188 !!$ IF( lendta(igrd) <= 0 ) THEN 189 !!$ idvar = iom_varid( inum,'u1' ) 190 !!$ lendta(igrd) = iom_file(inum)%dimsz(1,idvar) 191 !!$ WRITE(numout,*) 'Dim size for u1 is ',lendta(igrd) 192 !!$ ENDIF 193 !!$ ipi = lendta(igrd) 194 !!$ CALL iom_get( inum, jpdom_unknown, 'u1', zdta(1:ipi,1:ipj) ) 195 !!$ DO ib = 1, nblenrim(igrd) 196 !!$ u1(ib,itide) = zdta(nbmap(ib,igrd),1) 197 !!$ END DO 198 !!$ CALL iom_get( inum, jpdom_unknown, 'u2', zdta(1:ipi,1:ipj) ) 199 !!$ DO ib = 1, nblenrim(igrd) 200 !!$ u2(ib,itide) = zdta(nbmap(ib,igrd),1) 201 !!$ END DO 202 !!$ CALL iom_close( inum ) 203 !!$ ! 204 !!$ ! ! V fields 205 !!$ clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc' 206 !!$ if(lwp) write(numout,*) 'Reading data from file ', clfile 207 !!$ CALL iom_open( clfile, inum ) 208 !!$ igrd = 6 209 !!$ IF( lendta(igrd) <= 0 ) THEN 210 !!$ idvar = iom_varid( inum,'v1' ) 211 !!$ lendta(igrd) = iom_file(inum)%dimsz(1,idvar) 212 !!$ WRITE(numout,*) 'Dim size for v1 is ', lendta(igrd) 213 !!$ ENDIF 214 !!$ ipi = lendta(igrd) 215 !!$ CALL iom_get( inum, jpdom_unknown, 'v1', zdta(1:ipi,1:ipj) ) 216 !!$ DO ib = 1, nblenrim(igrd) 217 !!$ v1(ib,itide) = zdta(nbmap(ib,igrd),1) 218 !!$ END DO 219 !!$ CALL iom_get( inum, jpdom_unknown, 'v2', zdta(1:ipi,1:ipj) ) 220 !!$ DO ib=1, nblenrim(igrd) 221 !!$ v2(ib,itide) = zdta(nbmap(ib,igrd),1) 222 !!$ END DO 223 !!$ CALL iom_close( inum ) 224 !!$ ! 225 !!$ END DO ! end loop on tidal components 226 !!$ 227 !!$ IF( ln_tide_date ) THEN ! correct for date factors 228 !!$ 229 !!$!! used nmonth, nyear and nday from daymod.... 230 !!$ ! Calculate date corrects for 15 standard consituents 231 !!$ ! This is the initialisation step, so nday, nmonth, nyear are the 232 !!$ ! initial date/time of the integration. 233 !!$ print *, nday,nmonth,nyear 234 !!$ nyear = int(ndate0 / 10000 ) ! initial year 235 !!$ nmonth = int((ndate0 - nyear * 10000 ) / 100 ) ! initial month 236 !!$ nday = int(ndate0 - nyear * 10000 - nmonth * 100) 237 !!$ 238 !!$ CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu ) 239 !!$ 240 !!$ IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear 241 !!$ 242 !!$ DO itide = 1, ntide ! loop on tidal components 243 !!$ ! 244 !!$ IF( nindx(itide) /= 0 ) THEN 245 !!$!!gm use rpi and rad global variable 246 !!$ z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0 247 !!$ z_atde=z_ftc(nindx(itide))*cos(z_arg) 248 !!$ z_btde=z_ftc(nindx(itide))*sin(z_arg) 249 !!$ IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), tide_speed(itide), & 250 !!$ & z_ftc(nindx(itide)), z_vplu(nindx(itide)) 251 !!$ ELSE 252 !!$ z_atde = 1.0_wp 253 !!$ z_btde = 0.0_wp 254 !!$ ENDIF 255 !!$ ! ! elevation 256 !!$ igrd = 4 257 !!$ DO ib = 1, nblenrim(igrd) 258 !!$ z1t = z_atde * ssh1(ib,itide) + z_btde * ssh2(ib,itide) 259 !!$ z2t = z_atde * ssh2(ib,itide) - z_btde * ssh1(ib,itide) 260 !!$ ssh1(ib,itide) = z1t 261 !!$ ssh2(ib,itide) = z2t 262 !!$ END DO 263 !!$ ! ! u 264 !!$ igrd = 5 265 !!$ DO ib = 1, nblenrim(igrd) 266 !!$ z1t = z_atde * u1(ib,itide) + z_btde * u2(ib,itide) 267 !!$ z2t = z_atde * u2(ib,itide) - z_btde * u1(ib,itide) 268 !!$ u1(ib,itide) = z1t 269 !!$ u2(ib,itide) = z2t 270 !!$ END DO 271 !!$ ! ! v 272 !!$ igrd = 6 273 !!$ DO ib = 1, nblenrim(igrd) 274 !!$ z1t = z_atde * v1(ib,itide) + z_btde * v2(ib,itide) 275 !!$ z2t = z_atde * v2(ib,itide) - z_btde * v1(ib,itide) 276 !!$ v1(ib,itide) = z1t 277 !!$ v2(ib,itide) = z2t 278 !!$ END DO 279 !!$ ! 280 !!$ END DO ! end loop on tidal components 281 !!$ ! 282 !!$ ENDIF ! date correction 283 !!$ ! 284 !!$ END SUBROUTINE tide_data 285 !!$ 286 !!$ 287 !!$ SUBROUTINE tide_update ( kt, jit ) 288 !!$ !!---------------------------------------------------------------------- 289 !!$ !! *** SUBROUTINE tide_update *** 290 !!$ !! 291 !!$ !! ** Purpose : - Add tidal forcing to ssh_obc, ubt_obc and vbt_obc arrays. 292 !!$ !! 293 !!$ !!---------------------------------------------------------------------- 294 !!$ INTEGER, INTENT( in ) :: kt ! Main timestep counter 295 !!$!!gm doctor jit ==> kit 296 !!$ INTEGER, INTENT( in ) :: jit ! Barotropic timestep counter (for timesplitting option) 297 !!$ !! 298 !!$ INTEGER :: itide, igrd, ib ! dummy loop indices 299 !!$ REAL(wp) :: z_arg, z_sarg ! 300 !!$ REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost 301 !!$ !!---------------------------------------------------------------------- 302 !!$ 303 !!$ ! Note tide phase speeds are in deg/hour, so we need to convert the 304 !!$ ! elapsed time in seconds to hours by dividing by 3600.0 305 !!$ IF( jit == 0 ) THEN 306 !!$ z_arg = kt * rdt * rad / 3600.0 307 !!$ ELSE ! we are in a barotropic subcycle (for timesplitting option) 308 !!$! z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,lwp) ) * rad / 3600.0 309 !!$ z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,wp) ) * rad / 3600.0 310 !!$ ENDIF 311 !!$ 312 !!$ DO itide = 1, ntide 313 !!$ z_sarg = z_arg * tide_speed(itide) 314 !!$ z_cost(itide) = COS( z_sarg ) 315 !!$ z_sist(itide) = SIN( z_sarg ) 316 !!$ END DO 317 !!$ 318 !!$ ! summing of tidal constituents into BDY arrays 319 !!$ sshtide(:) = 0.0 320 !!$ utide (:) = 0.0 321 !!$ vtide (:) = 0.0 322 !!$ ! 323 !!$ DO itide = 1, ntide 324 !!$ igrd=4 ! SSH on tracer grid. 325 !!$ DO ib = 1, nblenrim(igrd) 326 !!$ sshtide(ib) =sshtide(ib)+ ssh1(ib,itide)*z_cost(itide) + ssh2(ib,itide)*z_sist(itide) 327 !!$ ! if(lwp) write(numout,*) 'z',ib,itide,sshtide(ib), ssh1(ib,itide),ssh2(ib,itide) 328 !!$ END DO 329 !!$ igrd=5 ! U grid 330 !!$ DO ib=1, nblenrim(igrd) 331 !!$ utide(ib) = utide(ib)+ u1(ib,itide)*z_cost(itide) + u2(ib,itide)*z_sist(itide) 332 !!$ ! if(lwp) write(numout,*) 'u',ib,itide,utide(ib), u1(ib,itide),u2(ib,itide) 333 !!$ END DO 334 !!$ igrd=6 ! V grid 335 !!$ DO ib=1, nblenrim(igrd) 336 !!$ vtide(ib) = vtide(ib)+ v1(ib,itide)*z_cost(itide) + v2(ib,itide)*z_sist(itide) 337 !!$ ! if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), v1(ib,itide),v2(ib,itide) 338 !!$ END DO 339 !!$ END DO 340 !!$ ! 341 !!$ END SUBROUTINE tide_update 342 !!$ 343 !!$!!gm doctor naming of dummy argument variables!!! and all local variables 344 !!$ 345 !!$ SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu ) 346 !!$ !!---------------------------------------------------------------------- 347 !!$ !! *** SUBROUTINE uvset *** 348 !!$ !! 349 !!$ !! ** Purpose : - adjust tidal forcing for date factors 350 !!$ !! 351 !!$ !!---------------------------------------------------------------------- 352 !!$ implicit none 353 !!$ INTEGER, INTENT( in ) :: ihs ! Start time hours 354 !!$ INTEGER, INTENT( in ) :: iday ! start time days 355 !!$ INTEGER, INTENT( in ) :: imnth ! start time month 356 !!$ INTEGER, INTENT( in ) :: iyr ! start time year 357 !!$ !! 358 !!$!!gm nc is jptides_max ???? 359 !!$ INTEGER , PARAMETER :: nc =15 ! maximum number of constituents 360 !!$ CHARACTER(len=8), DIMENSION(nc) :: cname 361 !!$ INTEGER :: year, vd, ivdy, ndc, i, k 362 !!$ REAL(wp) :: ss, h, p, en, p1, rtd 363 !!$ REAL(wp), DIMENSION(nc) :: f ! nodal correction 364 !!$ REAL(wp), DIMENSION(nc) :: z_vplu ! phase correction 365 !!$ REAL(wp), DIMENSION(nc) :: u, v, zig 366 !!$ !! 367 !!$ DATA cname/ 'q1' , 'o1' , 'p1' , 's1' , 'k1' , & 368 !!$ & '2n2' , 'mu2' , 'n2' , 'nu2' , 'm2' , & 369 !!$ & 'l2' , 't2' , 's2' , 'k2' , 'm4' / 370 !!$ DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878, .2625161701, & 371 !!$ & .4868657873, .4881373225, .4963669182, .4976384533, .5058680490, & 372 !!$ & .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098 / 373 !!$ !!---------------------------------------------------------------------- 374 !!$! 375 !!$! ihs - start time gmt on ... 376 !!$! iday/imnth/iyr - date e.g. 12/10/87 377 !!$! 378 !!$ CALL vday(iday,imnth,iyr,ivdy) 379 !!$ vd = ivdy 380 !!$! 381 !!$!rp note change of year number for d. blackman shpen 382 !!$!rp if(iyr.ge.1000) year=iyr-1900 383 !!$!rp if(iyr.lt.1000) year=iyr 384 !!$ year = iyr 385 !!$! 386 !!$!.....year = year of required data 387 !!$!.....vd = day of required data..set up for 0000gmt day year 388 !!$ ndc = nc 389 !!$!.....ndc = number of constituents allowed 390 !!$!!gm use rpi ? 391 !!$ rtd = 360.0 / 6.2831852 392 !!$ DO i = 1, ndc 393 !!$ zig(i) = zig(i)*rtd 394 !!$ ! sigo(i)= zig(i) 395 !!$ END DO 396 !!$ 397 !!$!!gm try to avoid RETURN in F90 398 !!$ IF( year == 0 ) RETURN 399 !!$ 400 !!$ CALL shpen( year, vd, ss, h , p , en, p1 ) 401 !!$ CALL ufset( p , en, u , f ) 402 !!$ CALL vset ( ss , h , p , en, p1, v ) 403 !!$ ! 404 !!$ DO k = 1, nc 405 !!$ z_vplu(k) = v(k) + u(k) 406 !!$ z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k) 407 !!$ DO WHILE( z_vplu(k) < 0 ) 408 !!$ z_vplu(k) = z_vplu(k) + 360.0 409 !!$ END DO 410 !!$ DO WHILE( z_vplu(k) > 360. ) 411 !!$ z_vplu(k) = z_vplu(k) - 360.0 412 !!$ END DO 413 !!$ END DO 414 !!$ ! 415 !!$ END SUBROUTINE uvset 416 !!$ 417 !!$ 418 !!$ SUBROUTINE vday( iday, imnth, iy, ivdy ) 419 !!$ !!---------------------------------------------------------------------- 420 !!$ !! *** SUBROUTINE vday *** 421 !!$ !! 422 !!$ !! ** Purpose : - adjust tidal forcing for date factors 423 !!$ !! 424 !!$ !!---------------------------------------------------------------------- 425 !!$ INTEGER, INTENT(in ) :: iday, imnth, iy ! ???? 426 !!$ INTEGER, INTENT( out) :: ivdy ! ??? 427 !!$ !! 428 !!$ INTEGER :: iyr 429 !!$ !!------------------------------------------------------------------------------ 430 !!$ 431 !!$!!gm nday_year in day mode is the variable compiuted here, no? 432 !!$!!gm nday_year , & !: curent day counted from jan 1st of the current year 433 !!$ 434 !!$ !calculate day number in year from day/month/year 435 !!$ if(imnth.eq.1) ivdy=iday 436 !!$ if(imnth.eq.2) ivdy=iday+31 437 !!$ if(imnth.eq.3) ivdy=iday+59 438 !!$ if(imnth.eq.4) ivdy=iday+90 439 !!$ if(imnth.eq.5) ivdy=iday+120 440 !!$ if(imnth.eq.6) ivdy=iday+151 441 !!$ if(imnth.eq.7) ivdy=iday+181 442 !!$ if(imnth.eq.8) ivdy=iday+212 443 !!$ if(imnth.eq.9) ivdy=iday+243 444 !!$ if(imnth.eq.10) ivdy=iday+273 445 !!$ if(imnth.eq.11) ivdy=iday+304 446 !!$ if(imnth.eq.12) ivdy=iday+334 447 !!$ iyr=iy 448 !!$ if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1 449 !!$ if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1 450 !!$ if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1 451 !!$ ! 452 !!$ END SUBROUTINE vday 453 !!$ 454 !!$!!doctor norme for dummy arguments 455 !!$ 456 !!$ SUBROUTINE shpen( year, vd, s, h, p, en, p1 ) 457 !!$ !!---------------------------------------------------------------------- 458 !!$ !! *** SUBROUTINE shpen *** 459 !!$ !! 460 !!$ !! ** Purpose : - calculate astronomical arguments for tides 461 !!$ !! this version from d. blackman 30 nove 1990 462 !!$ !! 463 !!$ !!---------------------------------------------------------------------- 464 !!$!!gm add INTENT in, out or inout.... DOCTOR name.... 465 !!$!!gm please do not use variable name with a single letter: impossible to search in a code 466 !!$ INTEGER :: year,vd 467 !!$ REAL(wp) :: s,h,p,en,p1 468 !!$ !! 469 !!$ INTEGER :: yr,ilc,icent,it,iday,ild,ipos,nn,iyd 470 !!$ REAL(wp) :: cycle,t,td,delt(84),delta,deltat 471 !!$ !! 472 !!$ DATA delt /-5.04, -3.90, -2.87, -0.58, 0.71, 1.80, & 473 !!$ & 3.08, 4.63, 5.86, 7.21, 8.58, 10.50, 12.10, & 474 !!$ & 12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39, & 475 !!$ & 19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68, & 476 !!$ & 22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63, & 477 !!$ & 23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80, & 478 !!$ & 24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89, & 479 !!$ & 27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96, & 480 !!$ & 31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39, & 481 !!$ & 33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87, & 482 !!$ & 38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00, & 483 !!$ & 45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81, & 484 !!$ & 52.57 / 485 !!$ !!---------------------------------------------------------------------- 486 !!$ 487 !!$ cycle = 360.0 488 !!$ ilc = 0 489 !!$ icent = year / 100 490 !!$ yr = year - icent * 100 491 !!$ t = icent - 20 492 !!$! 493 !!$! for the following equations 494 !!$! time origin is fixed at 00 hr of jan 1st,2000. 495 !!$! see notes by cartwright 496 !!$! 497 !!$!!gm old coding style, use CASE instead and avoid GOTO (obsolescence in fortran 90) 498 !!$!!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 499 !!$ it = icent - 20 500 !!$ if (it) 1,2,2 501 !!$ 1 iday = it/4 -it 502 !!$ go to 3 503 !!$ 2 iday = (it+3)/4 - it 504 !!$! 505 !!$! t is in julian century 506 !!$! correction in gegorian calander where only century year divisible 507 !!$! by 4 is leap year. 508 !!$! 509 !!$ 3 continue 510 !!$! 511 !!$ td = 0.0 512 !!$! 513 !!$!!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 514 !!$ if (yr) 4,5,4 515 !!$! 516 !!$ 4 iyd = 365*yr 517 !!$ ild = (yr-1)/4 518 !!$ if((icent - (icent/4)*4) .eq. 0) ilc = 1 519 !!$ td = iyd + ild + ilc 520 !!$! 521 !!$ 5 td = td + iday + vd -1.0 - 0.5 522 !!$ t = t + (td/36525.0) 523 !!$! 524 !!$ ipos=year-1899 525 !!$ if (ipos .lt. 0) go to 7 526 !!$ if (ipos .gt. 83) go to 6 527 !!$! 528 !!$ delta = (delt(ipos+1)+delt(ipos))/2.0 529 !!$ go to 7 530 !!$! 531 !!$ 6 delta= (65.0-50.5)/20.0*(year-1980)+50.5 532 !!$! 533 !!$ 7 deltat = delta * 1.0e-6 534 !!$! 535 !!$!!gm precision of the computation : example for s it should be replace by: 536 !!$!!gm s = 218.3165 + (481267.8813 - 0.0016*t)*t + 152.0*deltat ==> more precise modify the last digits results 537 !!$ s = 218.3165 + 481267.8813*t - 0.0016*t*t + 152.0*deltat 538 !!$ h = 280.4661 + 36000.7698 *t + 0.0003*t*t + 11.0*deltat 539 !!$ p = 83.3535 + 4069.0139 *t - 0.0103*t*t + deltat 540 !!$ en = 234.9555 + 1934.1363 *t - 0.0021*t*t + deltat 541 !!$ p1 = 282.9384 + 1.7195 *t + 0.0005*t*t 542 !!$ ! 543 !!$ nn = s / cycle 544 !!$ s = s - nn * cycle 545 !!$ IF( s < 0.e0 ) s = s + cycle 546 !!$ ! 547 !!$ nn = h / cycle 548 !!$ h = h - cycle * nn 549 !!$ IF( h < 0.e0 ) h = h + cycle 550 !!$ ! 551 !!$ nn = p / cycle 552 !!$ p = p - cycle * nn 553 !!$ IF( p < 0.e0) p = p + cycle 554 !!$ ! 555 !!$ nn = en / cycle 556 !!$ en = en - cycle * nn 557 !!$ IF( en < 0.e0 ) en = en + cycle 558 !!$ en = cycle - en 559 !!$ ! 560 !!$ nn = p1 / cycle 561 !!$ p1 = p1 - nn * cycle 562 !!$ ! 563 !!$ END SUBROUTINE shpen 564 !!$ 565 !!$ 566 !!$ SUBROUTINE ufset( p, cn, b, a ) 567 !!$ !!---------------------------------------------------------------------- 568 !!$ !! *** SUBROUTINE ufset *** 569 !!$ !! 570 !!$ !! ** Purpose : - calculate nodal parameters for the tides 571 !!$ !! 572 !!$ !!---------------------------------------------------------------------- 573 !!$!!gm doctor naming of dummy argument variables!!! and all local variables 574 !!$!!gm nc is jptides_max ???? 575 !!$ integer nc 576 !!$ parameter (nc=15) 577 !!$ REAL(wp) p,cn 578 !!$ !! 579 !!$ 580 !!$!!gm rad is already a public variable defined in phycst.F90 .... ==> doctor norme local real start with "z" 581 !!$ REAL(wp) :: w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad 582 !!$ REAL(wp) :: a(nc), b(nc) 583 !!$ INTEGER :: ndc, k 584 !!$ !!---------------------------------------------------------------------- 585 !!$ 586 !!$ ndc = nc 587 !!$ 588 !!$! a=f , b =u 589 !!$! t is zero as compared to tifa. 590 !!$!! use rad defined in phycst (i.e. add a USE phycst at the begining of the module 591 !!$ rad = 6.2831852d0/360.0 592 !!$ pw = p * rad 593 !!$ nw = cn * rad 594 !!$ w1 = cos( nw ) 595 !!$ w2 = cos( 2*nw ) 596 !!$ w3 = cos( 3*nw ) 597 !!$ w4 = sin( nw ) 598 !!$ w5 = sin( 2*nw ) 599 !!$ w6 = sin( 3*nw ) 600 !!$ w7 = 1. - 0.2505 * COS( 2*pw ) - 0.1102 * COS(2*pw-nw ) & 601 !!$ & - 0.156 * COS( 2*pw-2*nw ) - 0.037 * COS( nw ) 602 !!$ w8 = - 0.2505 * SIN( 2*pw ) - 0.1102 * SIN(2*pw-nw ) & 603 !!$ & - 0.0156 * SIN( 2*pw-2*nw ) - 0.037 * SIN( nw ) 604 !!$ ! 605 !!$ a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3 606 !!$ b(1) = 0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6 607 !!$ ! q1 608 !!$ a(2) = a(1) 609 !!$ b(2) = b(1) 610 !!$ ! o1 611 !!$ a(3) = 1.0 612 !!$ b(3) = 0.0 613 !!$ ! p1 614 !!$ a(4) = 1.0 615 !!$ b(4) = 0.0 616 !!$ ! s1 617 !!$ a(5) = 1.0060+0.1150*w1- 0.0088*w2 +0.0006*w3 618 !!$ b(5) = -0.1546*w4 + 0.0119*w5 -0.0012*w6 619 !!$ ! k1 620 !!$ a(6) =1.0004 -0.0373*w1+ 0.0002*w2 621 !!$ b(6) = -0.0374*w4 622 !!$ ! 2n2 623 !!$ a(7) = a(6) 624 !!$ b(7) = b(6) 625 !!$ ! mu2 626 !!$ a(8) = a(6) 627 !!$ b(8) = b(6) 628 !!$ ! n2 629 !!$ a(9) = a(6) 630 !!$ b(9) = b(6) 631 !!$ ! nu2 632 !!$ a(10) = a(6) 633 !!$ b(10) = b(6) 634 !!$ ! m2 635 !!$ a(11) = SQRT( w7 * w7 + w8 * w8 ) 636 !!$ b(11) = ATAN( w8 / w7 ) 637 !!$!!gmuse rpi instead of 3.141992 ??? true pi is rpi=3.141592653589793_wp ..... ???? 638 !!$ IF( w7 < 0.e0 ) b(11) = b(11) + 3.141992 639 !!$ ! l2 640 !!$ a(12) = 1.0 641 !!$ b(12) = 0.0 642 !!$ ! t2 643 !!$ a(13)= a(12) 644 !!$ b(13)= b(12) 645 !!$ ! s2 646 !!$ a(14) = 1.0241+0.2863*w1+0.0083*w2 -0.0015*w3 647 !!$ b(14) = -0.3096*w4 + 0.0119*w5 - 0.0007*w6 648 !!$ ! k2 649 !!$ a(15) = a(6)*a(6) 650 !!$ b(15) = 2*b(6) 651 !!$ ! m4 652 !!$!!gm old coding, remove GOTO and label of lines 653 !!$!!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 654 !!$ DO 40 k = 1,ndc 655 !!$ b(k) = b(k)/rad 656 !!$32 if (b(k)) 34,35,35 657 !!$34 b(k) = b(k) + 360.0 658 !!$ go to 32 659 !!$35 if (b(k)-360.0) 40,37,37 660 !!$37 b(k) = b(k)-360.0 661 !!$ go to 35 662 !!$40 continue 663 !!$ ! 664 !!$ END SUBROUTINE ufset 665 !!$ 666 !!$ 667 !!$ SUBROUTINE vset( s,h,p,en,p1,v) 668 !!$ !!---------------------------------------------------------------------- 669 !!$ !! *** SUBROUTINE vset *** 670 !!$ !! 671 !!$ !! ** Purpose : - calculate tidal phases for 0000gmt on start day of run 672 !!$ !! 673 !!$ !!---------------------------------------------------------------------- 674 !!$!!gm doctor naming of dummy argument variables!!! and all local variables 675 !!$!!gm nc is jptides_max ???? 676 !!$!!gm en argument is not used: suppress it ? 677 !!$ integer nc 678 !!$ parameter (nc=15) 679 !!$ real(wp) s,h,p,en,p1 680 !!$ real(wp) v(nc) 681 !!$ !! 682 !!$ integer ndc, k 683 !!$ !!---------------------------------------------------------------------- 684 !!$ 685 !!$ ndc = nc 686 !!$ ! v s are computed here. 687 !!$ v(1) =-3*s +h +p +270 ! Q1 688 !!$ v(2) =-2*s +h +270.0 ! O1 689 !!$ v(3) =-h +270 ! P1 690 !!$ v(4) =180 ! S1 691 !!$ v(5) =h +90.0 ! K1 692 !!$ v(6) =-4*s +2*h +2*p ! 2N2 693 !!$ v(7) =-4*(s-h) ! MU2 694 !!$ v(8) =-3*s +2*h +p ! N2 695 !!$ v(9) =-3*s +4*h -p ! MU2 696 !!$ v(10) =-2*s +2*h ! M2 697 !!$ v(11) =-s +2*h -p +180 ! L2 698 !!$ v(12) =-h +p1 ! T2 699 !!$ v(13) =0 ! S2 700 !!$ v(14) =h+h ! K2 701 !!$ v(15) =2*v(10) ! M4 702 !!$! 703 !!$!!gm old coding, remove GOTO and label of lines 704 !!$!!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 705 !!$ do 72 k = 1, ndc 706 !!$69 if( v(k) ) 70,71,71 707 !!$70 v(k) = v(k)+360.0 708 !!$ go to 69 709 !!$71 if( v(k) - 360.0 ) 72,73,73 710 !!$73 v(k) = v(k)-360.0 711 !!$ go to 71 712 !!$72 continue 713 !!$ ! 714 !!$ END SUBROUTINE vset 715 !!$ 716 !!$#else 12 !!---------------------------------------------------------------------- 13 !! 'key_obc' Open Boundary Condition 14 !!---------------------------------------------------------------------- 15 !! PUBLIC 16 !! tide_init : read of namelist and initialisation of tidal harmonics data 17 !! tide_update : calculation of tidal forcing at each timestep 18 !! PRIVATE 19 !! uvset :\ 20 !! vday : | Routines to correct tidal harmonics forcing for 21 !! shpen : | start time of integration 22 !! ufset : | 23 !! vset :/ 24 !!---------------------------------------------------------------------- 25 USE oce ! ocean dynamics and tracers 26 USE dom_oce ! ocean space and time domain 27 USE iom 28 USE in_out_manager ! I/O units 29 USE phycst ! physical constants 30 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 31 USE obc_par ! Unstructured boundary parameters 32 USE obc_oce ! ocean open boundary conditions 33 USE daymod ! calendar 34 USE fldread, ONLY: fld_map 35 36 IMPLICIT NONE 37 PRIVATE 38 39 PUBLIC tide_init ! routine called in nemo_init 40 PUBLIC tide_update ! routine called in obcdyn 41 42 TYPE, PUBLIC :: TIDES_DATA !: Storage for external tidal harmonics data 43 INTEGER :: ncpt !: Actual number of tidal components 44 REAL(wp), POINTER, DIMENSION(:) :: speed !: Phase speed of tidal constituent (deg/hr) 45 REAL(wp), POINTER, DIMENSION(:,:,:) :: ssh !: Tidal constituents : SSH 46 REAL(wp), POINTER, DIMENSION(:,:,:) :: u !: Tidal constituents : U 47 REAL(wp), POINTER, DIMENSION(:,:,:) :: v !: Tidal constituents : V 48 END TYPE TIDES_DATA 49 50 INTEGER, PUBLIC, PARAMETER :: jptides_max = 15 !: Max number of tidal contituents 51 52 TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_obc), TARGET :: tides !: External tidal harmonics data 53 54 !!---------------------------------------------------------------------- 55 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 56 !! $Id: obctides.F90 2528 2010-12-27 17:33:53Z rblod $ 57 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 58 !!---------------------------------------------------------------------- 59 CONTAINS 60 61 SUBROUTINE tide_init 62 !!---------------------------------------------------------------------- 63 !! *** SUBROUTINE tide_init *** 64 !! 65 !! ** Purpose : - Read in namelist for tides and initialise external 66 !! tidal harmonics data 67 !! 68 !!---------------------------------------------------------------------- 69 !! namelist variables 70 !!------------------- 71 LOGICAL :: ln_tide_date !: =T correct tide phases and amplitude for model start date 72 CHARACTER(len=80) :: filtide !: Filename root for tidal input files 73 CHARACTER(len= 4), DIMENSION(jptides_max) :: tide_cpt !: Names of tidal components used. 74 REAL(wp), DIMENSION(jptides_max) :: tide_speed !: Phase speed of tidal constituent (deg/hr) 75 !! 76 INTEGER, DIMENSION(jptides_max) :: nindx !: index of pre-set tidal components 77 INTEGER :: ib_obc, itide, ib !: dummy loop indices 78 INTEGER :: inum, igrd 79 INTEGER, DIMENSION(3) :: ilen0 !: length of boundary data (from OBC arrays) 80 CHARACTER(len=80) :: clfile !: full file name for tidal input file 81 REAL(wp) :: z_arg, z_atde, z_btde, z1t, z2t 82 REAL(wp),DIMENSION(jptides_max) :: z_vplu, z_ftc 83 REAL(wp),ALLOCATABLE, DIMENSION(:,:,:) :: dta_read !: work space to read in tidal harmonics data 84 !! 85 TYPE(TIDES_DATA), POINTER :: td !: local short cut 86 !! 87 NAMELIST/namobc_tide/ln_tide_date, filtide, tide_cpt, tide_speed 88 !!---------------------------------------------------------------------- 89 90 IF(lwp) WRITE(numout,*) 91 IF(lwp) WRITE(numout,*) 'tide_init : initialization of tidal harmonic forcing at open boundaries' 92 IF(lwp) WRITE(numout,*) '~~~~~~~~~' 93 94 REWIND(numnam) 95 DO ib_obc = 1, nb_obc 96 IF( nn_tides(ib_obc) .gt. 0 ) THEN 97 98 td => tides(ib_obc) 99 100 ! Namelist namobc_tide : tidal harmonic forcing at open boundaries 101 ln_tide_date = .false. 102 filtide(:) = '' 103 tide_speed(:) = 0.0 104 tide_cpt(:) = '' 105 106 ! Don't REWIND here - may need to read more than one of these namelists. 107 READ ( numnam, namobc_tide ) 108 ! ! Count number of components specified 109 td%ncpt = 0 110 DO itide = 1, jptides_max 111 IF( tide_cpt(itide) /= '' ) THEN 112 td%ncpt = td%ncpt + 1 113 ENDIF 114 END DO 115 116 ! Fill in phase speeds from namelist 117 ALLOCATE( td%speed(td%ncpt) ) 118 td%speed = tide_speed(1:td%ncpt) 119 120 ! Find constituents in standard list 121 DO itide = 1, td%ncpt 122 nindx(itide) = 0 123 IF( TRIM( tide_cpt(itide) ) == 'Q1' ) nindx(itide) = 1 124 IF( TRIM( tide_cpt(itide) ) == 'O1' ) nindx(itide) = 2 125 IF( TRIM( tide_cpt(itide) ) == 'P1' ) nindx(itide) = 3 126 IF( TRIM( tide_cpt(itide) ) == 'S1' ) nindx(itide) = 4 127 IF( TRIM( tide_cpt(itide) ) == 'K1' ) nindx(itide) = 5 128 IF( TRIM( tide_cpt(itide) ) == '2N2' ) nindx(itide) = 6 129 IF( TRIM( tide_cpt(itide) ) == 'MU2' ) nindx(itide) = 7 130 IF( TRIM( tide_cpt(itide) ) == 'N2' ) nindx(itide) = 8 131 IF( TRIM( tide_cpt(itide) ) == 'NU2' ) nindx(itide) = 9 132 IF( TRIM( tide_cpt(itide) ) == 'M2' ) nindx(itide) = 10 133 IF( TRIM( tide_cpt(itide) ) == 'L2' ) nindx(itide) = 11 134 IF( TRIM( tide_cpt(itide) ) == 'T2' ) nindx(itide) = 12 135 IF( TRIM( tide_cpt(itide) ) == 'S2' ) nindx(itide) = 13 136 IF( TRIM( tide_cpt(itide) ) == 'K2' ) nindx(itide) = 14 137 IF( TRIM( tide_cpt(itide) ) == 'M4' ) nindx(itide) = 15 138 IF( nindx(itide) == 0 .AND. lwp ) THEN 139 WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list' 140 CALL ctl_warn( ctmp1 ) 141 ENDIF 142 END DO 143 ! ! Parameter control and print 144 IF( td%ncpt < 1 ) THEN 145 CALL ctl_stop( ' Did not find any tidal components in namelist namobc_tide' ) 146 ELSE 147 IF(lwp) WRITE(numout,*) ' Namelist namobc_tide : tidal harmonic forcing at open boundaries' 148 IF(lwp) WRITE(numout,*) ' tidal components specified ', td%ncpt 149 IF(lwp) WRITE(numout,*) ' ', tide_cpt(1:td%ncpt) 150 IF(lwp) WRITE(numout,*) ' associated phase speeds (deg/hr) : ' 151 IF(lwp) WRITE(numout,*) ' ', tide_speed(1:td%ncpt) 152 ENDIF 153 154 155 ! Allocate space for tidal harmonics data - 156 ! get size from OBC data arrays 157 ! --------------------------------------- 158 159 ilen0(1) = SIZE( dta_obc(ib_obc)%ssh ) 160 ALLOCATE( td%ssh( ilen0(1), td%ncpt, 2 ) ) 161 162 ilen0(2) = SIZE( dta_obc(ib_obc)%u2d ) 163 ALLOCATE( td%u( ilen0(2), td%ncpt, 2 ) ) 164 165 ilen0(3) = SIZE( dta_obc(ib_obc)%v2d ) 166 ALLOCATE( td%v( ilen0(3), td%ncpt, 2 ) ) 167 168 ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) ) 169 170 171 ! Open files and read in tidal forcing data 172 ! ----------------------------------------- 173 174 DO itide = 1, td%ncpt 175 ! ! SSH fields 176 clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc' 177 IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 178 CALL iom_open( clfile, inum ) 179 CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,1) ) 180 td%ssh(:,itide,1) = dta_read(1:ilen0(1),1,1) 181 CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,1) ) 182 td%ssh(:,itide,2) = dta_read(1:ilen0(1),1,1) 183 CALL iom_close( inum ) 184 ! ! U fields 185 clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc' 186 IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 187 CALL iom_open( clfile, inum ) 188 CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,2) ) 189 td%u(:,itide,1) = dta_read(1:ilen0(2),1,1) 190 CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,2) ) 191 td%u(:,itide,2) = dta_read(1:ilen0(2),1,1) 192 CALL iom_close( inum ) 193 ! ! V fields 194 clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc' 195 IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 196 CALL iom_open( clfile, inum ) 197 CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,3) ) 198 td%v(:,itide,1) = dta_read(1:ilen0(3),1,1) 199 CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,3) ) 200 td%v(:,itide,2) = dta_read(1:ilen0(3),1,1) 201 CALL iom_close( inum ) 202 ! 203 END DO ! end loop on tidal components 204 205 IF( ln_tide_date ) THEN ! correct for date factors 206 207 !! used nmonth, nyear and nday from daymod.... 208 ! Calculate date corrects for 15 standard consituents 209 ! This is the initialisation step, so nday, nmonth, nyear are the 210 ! initial date/time of the integration. 211 print *, nday,nmonth,nyear 212 nyear = int(ndate0 / 10000 ) ! initial year 213 nmonth = int((ndate0 - nyear * 10000 ) / 100 ) ! initial month 214 nday = int(ndate0 - nyear * 10000 - nmonth * 100) 215 216 CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu ) 217 218 IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear 219 220 DO itide = 1, td%ncpt ! loop on tidal components 221 ! 222 IF( nindx(itide) /= 0 ) THEN 223 !!gm use rpi and rad global variable 224 z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0 225 z_atde=z_ftc(nindx(itide))*cos(z_arg) 226 z_btde=z_ftc(nindx(itide))*sin(z_arg) 227 IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), td%speed(itide), & 228 & z_ftc(nindx(itide)), z_vplu(nindx(itide)) 229 ELSE 230 z_atde = 1.0_wp 231 z_btde = 0.0_wp 232 ENDIF 233 ! ! elevation 234 igrd = 1 235 DO ib = 1, ilen0(igrd) 236 z1t = z_atde * td%ssh(ib,itide,1) + z_btde * td%ssh(ib,itide,2) 237 z2t = z_atde * td%ssh(ib,itide,2) - z_btde * td%ssh(ib,itide,1) 238 td%ssh(ib,itide,1) = z1t 239 td%ssh(ib,itide,2) = z2t 240 END DO 241 ! ! u 242 igrd = 2 243 DO ib = 1, ilen0(igrd) 244 z1t = z_atde * td%u(ib,itide,1) + z_btde * td%u(ib,itide,2) 245 z2t = z_atde * td%u(ib,itide,2) - z_btde * td%u(ib,itide,1) 246 td%u(ib,itide,1) = z1t 247 td%u(ib,itide,2) = z2t 248 END DO 249 ! ! v 250 igrd = 3 251 DO ib = 1, ilen0(igrd) 252 z1t = z_atde * td%v(ib,itide,1) + z_btde * td%v(ib,itide,2) 253 z2t = z_atde * td%v(ib,itide,2) - z_btde * td%v(ib,itide,1) 254 td%v(ib,itide,1) = z1t 255 td%v(ib,itide,2) = z2t 256 END DO 257 ! 258 END DO ! end loop on tidal components 259 ! 260 ENDIF ! date correction 261 ! 262 ENDIF ! nn_tides .gt. 0 263 ! 264 END DO ! loop on ib_obc 265 266 END SUBROUTINE tide_init 267 268 269 SUBROUTINE tide_update ( kt, jit, idx, dta, td ) 270 !!---------------------------------------------------------------------- 271 !! *** SUBROUTINE tide_update *** 272 !! 273 !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays. 274 !! 275 !!---------------------------------------------------------------------- 276 INTEGER, INTENT( in ) :: kt ! Main timestep counter 277 !!gm doctor jit ==> kit 278 INTEGER, INTENT( in ) :: jit ! Barotropic timestep counter (for timesplitting option) 279 TYPE(OBC_INDEX), INTENT( in ) :: idx ! OBC indices 280 TYPE(OBC_DATA), INTENT(inout) :: dta ! OBC external data 281 TYPE(TIDES_DATA),INTENT( in ) :: td ! tidal harmonics data 282 !! 283 INTEGER :: itide, igrd, ib ! dummy loop indices 284 REAL(wp) :: z_arg, z_sarg 285 REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost 286 !!---------------------------------------------------------------------- 287 288 ! Note tide phase speeds are in deg/hour, so we need to convert the 289 ! elapsed time in seconds to hours by dividing by 3600.0 290 IF( jit == 0 ) THEN 291 z_arg = kt * rdt * rad / 3600.0 292 ELSE ! we are in a barotropic subcycle (for timesplitting option) 293 ! z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,lwp) ) * rad / 3600.0 294 z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,wp) ) * rad / 3600.0 295 ENDIF 296 297 DO itide = 1, td%ncpt 298 z_sarg = z_arg * td%speed(itide) 299 z_cost(itide) = COS( z_sarg ) 300 z_sist(itide) = SIN( z_sarg ) 301 END DO 302 303 ! 304 DO itide = 1, td%ncpt 305 igrd=1 ! SSH on tracer grid. 306 DO ib = 1, idx%nblenrim(igrd) 307 dta%ssh(ib) = dta%ssh(ib) + td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide) 308 ! if(lwp) write(numout,*) 'z', ib, itide, dta%ssh(ib), td%ssh(ib,itide,1),td%ssh(ib,itide,2) 309 END DO 310 igrd=2 ! U grid 311 DO ib=1, idx%nblenrim(igrd) 312 dta%u2d(ib) = dta%u2d(ib) + td%u(ib,itide,1)*z_cost(itide) + td%u(ib,itide,2)*z_sist(itide) 313 ! if(lwp) write(numout,*) 'u',ib,itide,utide(ib), td%u(ib,itide,1),td%u(ib,itide,2) 314 END DO 315 igrd=3 ! V grid 316 DO ib=1, idx%nblenrim(igrd) 317 dta%v2d(ib) = dta%v2d(ib) + td%v(ib,itide,1)*z_cost(itide) + td%v(ib,itide,2)*z_sist(itide) 318 ! if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), td%v(ib,itide,1),td%v(ib,itide,2) 319 END DO 320 END DO 321 ! 322 END SUBROUTINE tide_update 323 324 !!gm doctor naming of dummy argument variables!!! and all local variables 325 326 SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu ) 327 !!---------------------------------------------------------------------- 328 !! *** SUBROUTINE uvset *** 329 !! 330 !! ** Purpose : - adjust tidal forcing for date factors 331 !! 332 !!---------------------------------------------------------------------- 333 implicit none 334 INTEGER, INTENT( in ) :: ihs ! Start time hours 335 INTEGER, INTENT( in ) :: iday ! start time days 336 INTEGER, INTENT( in ) :: imnth ! start time month 337 INTEGER, INTENT( in ) :: iyr ! start time year 338 !! 339 !!gm nc is jptides_max ???? 340 INTEGER , PARAMETER :: nc =15 ! maximum number of constituents 341 CHARACTER(len=8), DIMENSION(nc) :: cname 342 INTEGER :: year, vd, ivdy, ndc, i, k 343 REAL(wp) :: ss, h, p, en, p1, rtd 344 REAL(wp), DIMENSION(nc) :: f ! nodal correction 345 REAL(wp), DIMENSION(nc) :: z_vplu ! phase correction 346 REAL(wp), DIMENSION(nc) :: u, v, zig 347 !! 348 DATA cname/ 'q1' , 'o1' , 'p1' , 's1' , 'k1' , & 349 & '2n2' , 'mu2' , 'n2' , 'nu2' , 'm2' , & 350 & 'l2' , 't2' , 's2' , 'k2' , 'm4' / 351 DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878, .2625161701, & 352 & .4868657873, .4881373225, .4963669182, .4976384533, .5058680490, & 353 & .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098 / 354 !!---------------------------------------------------------------------- 355 ! 356 ! ihs - start time gmt on ... 357 ! iday/imnth/iyr - date e.g. 12/10/87 358 ! 359 CALL vday(iday,imnth,iyr,ivdy) 360 vd = ivdy 361 ! 362 !rp note change of year number for d. blackman shpen 363 !rp if(iyr.ge.1000) year=iyr-1900 364 !rp if(iyr.lt.1000) year=iyr 365 year = iyr 366 ! 367 !.....year = year of required data 368 !.....vd = day of required data..set up for 0000gmt day year 369 ndc = nc 370 !.....ndc = number of constituents allowed 371 !!gm use rpi ? 372 rtd = 360.0 / 6.2831852 373 DO i = 1, ndc 374 zig(i) = zig(i)*rtd 375 ! sigo(i)= zig(i) 376 END DO 377 378 !!gm try to avoid RETURN in F90 379 IF( year == 0 ) RETURN 380 381 CALL shpen( year, vd, ss, h , p , en, p1 ) 382 CALL ufset( p , en, u , f ) 383 CALL vset ( ss , h , p , en, p1, v ) 384 ! 385 DO k = 1, nc 386 z_vplu(k) = v(k) + u(k) 387 z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k) 388 DO WHILE( z_vplu(k) < 0 ) 389 z_vplu(k) = z_vplu(k) + 360.0 390 END DO 391 DO WHILE( z_vplu(k) > 360. ) 392 z_vplu(k) = z_vplu(k) - 360.0 393 END DO 394 END DO 395 ! 396 END SUBROUTINE uvset 397 398 399 SUBROUTINE vday( iday, imnth, iy, ivdy ) 400 !!---------------------------------------------------------------------- 401 !! *** SUBROUTINE vday *** 402 !! 403 !! ** Purpose : - adjust tidal forcing for date factors 404 !! 405 !!---------------------------------------------------------------------- 406 INTEGER, INTENT(in ) :: iday, imnth, iy ! ???? 407 INTEGER, INTENT( out) :: ivdy ! ??? 408 !! 409 INTEGER :: iyr 410 !!------------------------------------------------------------------------------ 411 412 !!gm nday_year in day mode is the variable compiuted here, no? 413 !!gm nday_year , & !: curent day counted from jan 1st of the current year 414 415 !calculate day number in year from day/month/year 416 if(imnth.eq.1) ivdy=iday 417 if(imnth.eq.2) ivdy=iday+31 418 if(imnth.eq.3) ivdy=iday+59 419 if(imnth.eq.4) ivdy=iday+90 420 if(imnth.eq.5) ivdy=iday+120 421 if(imnth.eq.6) ivdy=iday+151 422 if(imnth.eq.7) ivdy=iday+181 423 if(imnth.eq.8) ivdy=iday+212 424 if(imnth.eq.9) ivdy=iday+243 425 if(imnth.eq.10) ivdy=iday+273 426 if(imnth.eq.11) ivdy=iday+304 427 if(imnth.eq.12) ivdy=iday+334 428 iyr=iy 429 if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1 430 if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1 431 if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1 432 ! 433 END SUBROUTINE vday 434 435 !!doctor norme for dummy arguments 436 437 SUBROUTINE shpen( year, vd, s, h, p, en, p1 ) 438 !!---------------------------------------------------------------------- 439 !! *** SUBROUTINE shpen *** 440 !! 441 !! ** Purpose : - calculate astronomical arguments for tides 442 !! this version from d. blackman 30 nove 1990 443 !! 444 !!---------------------------------------------------------------------- 445 !!gm add INTENT in, out or inout.... DOCTOR name.... 446 !!gm please do not use variable name with a single letter: impossible to search in a code 447 INTEGER :: year,vd 448 REAL(wp) :: s,h,p,en,p1 449 !! 450 INTEGER :: yr,ilc,icent,it,iday,ild,ipos,nn,iyd 451 REAL(wp) :: cycle,t,td,delt(84),delta,deltat 452 !! 453 DATA delt /-5.04, -3.90, -2.87, -0.58, 0.71, 1.80, & 454 & 3.08, 4.63, 5.86, 7.21, 8.58, 10.50, 12.10, & 455 & 12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39, & 456 & 19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68, & 457 & 22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63, & 458 & 23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80, & 459 & 24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89, & 460 & 27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96, & 461 & 31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39, & 462 & 33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87, & 463 & 38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00, & 464 & 45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81, & 465 & 52.57 / 466 !!---------------------------------------------------------------------- 467 468 cycle = 360.0 469 ilc = 0 470 icent = year / 100 471 yr = year - icent * 100 472 t = icent - 20 473 ! 474 ! for the following equations 475 ! time origin is fixed at 00 hr of jan 1st,2000. 476 ! see notes by cartwright 477 ! 478 !!gm old coding style, use CASE instead and avoid GOTO (obsolescence in fortran 90) 479 !!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 480 it = icent - 20 481 if (it) 1,2,2 482 1 iday = it/4 -it 483 go to 3 484 2 iday = (it+3)/4 - it 485 ! 486 ! t is in julian century 487 ! correction in gegorian calander where only century year divisible 488 ! by 4 is leap year. 489 ! 490 3 continue 491 ! 492 td = 0.0 493 ! 494 !!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 495 if (yr) 4,5,4 496 ! 497 4 iyd = 365*yr 498 ild = (yr-1)/4 499 if((icent - (icent/4)*4) .eq. 0) ilc = 1 500 td = iyd + ild + ilc 501 ! 502 5 td = td + iday + vd -1.0 - 0.5 503 t = t + (td/36525.0) 504 ! 505 ipos=year-1899 506 if (ipos .lt. 0) go to 7 507 if (ipos .gt. 83) go to 6 508 ! 509 delta = (delt(ipos+1)+delt(ipos))/2.0 510 go to 7 511 ! 512 6 delta= (65.0-50.5)/20.0*(year-1980)+50.5 513 ! 514 7 deltat = delta * 1.0e-6 515 ! 516 !!gm precision of the computation : example for s it should be replace by: 517 !!gm s = 218.3165 + (481267.8813 - 0.0016*t)*t + 152.0*deltat ==> more precise modify the last digits results 518 s = 218.3165 + 481267.8813*t - 0.0016*t*t + 152.0*deltat 519 h = 280.4661 + 36000.7698 *t + 0.0003*t*t + 11.0*deltat 520 p = 83.3535 + 4069.0139 *t - 0.0103*t*t + deltat 521 en = 234.9555 + 1934.1363 *t - 0.0021*t*t + deltat 522 p1 = 282.9384 + 1.7195 *t + 0.0005*t*t 523 ! 524 nn = s / cycle 525 s = s - nn * cycle 526 IF( s < 0.e0 ) s = s + cycle 527 ! 528 nn = h / cycle 529 h = h - cycle * nn 530 IF( h < 0.e0 ) h = h + cycle 531 ! 532 nn = p / cycle 533 p = p - cycle * nn 534 IF( p < 0.e0) p = p + cycle 535 ! 536 nn = en / cycle 537 en = en - cycle * nn 538 IF( en < 0.e0 ) en = en + cycle 539 en = cycle - en 540 ! 541 nn = p1 / cycle 542 p1 = p1 - nn * cycle 543 ! 544 END SUBROUTINE shpen 545 546 547 SUBROUTINE ufset( p, cn, b, a ) 548 !!---------------------------------------------------------------------- 549 !! *** SUBROUTINE ufset *** 550 !! 551 !! ** Purpose : - calculate nodal parameters for the tides 552 !! 553 !!---------------------------------------------------------------------- 554 !!gm doctor naming of dummy argument variables!!! and all local variables 555 !!gm nc is jptides_max ???? 556 integer nc 557 parameter (nc=15) 558 REAL(wp) p,cn 559 !! 560 561 !!gm rad is already a public variable defined in phycst.F90 .... ==> doctor norme local real start with "z" 562 REAL(wp) :: w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad 563 REAL(wp) :: a(nc), b(nc) 564 INTEGER :: ndc, k 565 !!---------------------------------------------------------------------- 566 567 ndc = nc 568 569 ! a=f , b =u 570 ! t is zero as compared to tifa. 571 !! use rad defined in phycst (i.e. add a USE phycst at the begining of the module 572 rad = 6.2831852d0/360.0 573 pw = p * rad 574 nw = cn * rad 575 w1 = cos( nw ) 576 w2 = cos( 2*nw ) 577 w3 = cos( 3*nw ) 578 w4 = sin( nw ) 579 w5 = sin( 2*nw ) 580 w6 = sin( 3*nw ) 581 w7 = 1. - 0.2505 * COS( 2*pw ) - 0.1102 * COS(2*pw-nw ) & 582 & - 0.156 * COS( 2*pw-2*nw ) - 0.037 * COS( nw ) 583 w8 = - 0.2505 * SIN( 2*pw ) - 0.1102 * SIN(2*pw-nw ) & 584 & - 0.0156 * SIN( 2*pw-2*nw ) - 0.037 * SIN( nw ) 585 ! 586 a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3 587 b(1) = 0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6 588 ! q1 589 a(2) = a(1) 590 b(2) = b(1) 591 ! o1 592 a(3) = 1.0 593 b(3) = 0.0 594 ! p1 595 a(4) = 1.0 596 b(4) = 0.0 597 ! s1 598 a(5) = 1.0060+0.1150*w1- 0.0088*w2 +0.0006*w3 599 b(5) = -0.1546*w4 + 0.0119*w5 -0.0012*w6 600 ! k1 601 a(6) =1.0004 -0.0373*w1+ 0.0002*w2 602 b(6) = -0.0374*w4 603 ! 2n2 604 a(7) = a(6) 605 b(7) = b(6) 606 ! mu2 607 a(8) = a(6) 608 b(8) = b(6) 609 ! n2 610 a(9) = a(6) 611 b(9) = b(6) 612 ! nu2 613 a(10) = a(6) 614 b(10) = b(6) 615 ! m2 616 a(11) = SQRT( w7 * w7 + w8 * w8 ) 617 b(11) = ATAN( w8 / w7 ) 618 !!gmuse rpi instead of 3.141992 ??? true pi is rpi=3.141592653589793_wp ..... ???? 619 IF( w7 < 0.e0 ) b(11) = b(11) + 3.141992 620 ! l2 621 a(12) = 1.0 622 b(12) = 0.0 623 ! t2 624 a(13)= a(12) 625 b(13)= b(12) 626 ! s2 627 a(14) = 1.0241+0.2863*w1+0.0083*w2 -0.0015*w3 628 b(14) = -0.3096*w4 + 0.0119*w5 - 0.0007*w6 629 ! k2 630 a(15) = a(6)*a(6) 631 b(15) = 2*b(6) 632 ! m4 633 !!gm old coding, remove GOTO and label of lines 634 !!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 635 DO 40 k = 1,ndc 636 b(k) = b(k)/rad 637 32 if (b(k)) 34,35,35 638 34 b(k) = b(k) + 360.0 639 go to 32 640 35 if (b(k)-360.0) 40,37,37 641 37 b(k) = b(k)-360.0 642 go to 35 643 40 continue 644 ! 645 END SUBROUTINE ufset 646 647 648 SUBROUTINE vset( s,h,p,en,p1,v) 649 !!---------------------------------------------------------------------- 650 !! *** SUBROUTINE vset *** 651 !! 652 !! ** Purpose : - calculate tidal phases for 0000gmt on start day of run 653 !! 654 !!---------------------------------------------------------------------- 655 !!gm doctor naming of dummy argument variables!!! and all local variables 656 !!gm nc is jptides_max ???? 657 !!gm en argument is not used: suppress it ? 658 integer nc 659 parameter (nc=15) 660 real(wp) s,h,p,en,p1 661 real(wp) v(nc) 662 !! 663 integer ndc, k 664 !!---------------------------------------------------------------------- 665 666 ndc = nc 667 ! v s are computed here. 668 v(1) =-3*s +h +p +270 ! Q1 669 v(2) =-2*s +h +270.0 ! O1 670 v(3) =-h +270 ! P1 671 v(4) =180 ! S1 672 v(5) =h +90.0 ! K1 673 v(6) =-4*s +2*h +2*p ! 2N2 674 v(7) =-4*(s-h) ! MU2 675 v(8) =-3*s +2*h +p ! N2 676 v(9) =-3*s +4*h -p ! MU2 677 v(10) =-2*s +2*h ! M2 678 v(11) =-s +2*h -p +180 ! L2 679 v(12) =-h +p1 ! T2 680 v(13) =0 ! S2 681 v(14) =h+h ! K2 682 v(15) =2*v(10) ! M4 683 ! 684 !!gm old coding, remove GOTO and label of lines 685 !!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 686 do 72 k = 1, ndc 687 69 if( v(k) ) 70,71,71 688 70 v(k) = v(k)+360.0 689 go to 69 690 71 if( v(k) - 360.0 ) 72,73,73 691 73 v(k) = v(k)-360.0 692 go to 71 693 72 continue 694 ! 695 END SUBROUTINE vset 696 697 #else 717 698 !!---------------------------------------------------------------------- 718 699 !! Dummy module NO Unstruct Open Boundary Conditions for tides -
branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
r2800 r2814 24 24 IMPLICIT NONE 25 25 PRIVATE 26 27 PUBLIC fld_map ! routine called by tides_init 26 28 27 29 TYPE, PUBLIC :: FLD_N !: Namelist field informations … … 629 631 !! using a general mapping (for open boundaries) 630 632 !!---------------------------------------------------------------------- 633 #if defined key_obc 631 634 USE obc_oce, ONLY: dta_global ! workspace to read in global data arrays 635 #endif 632 636 633 637 INTEGER , INTENT(in ) :: num ! stream number 634 638 CHARACTER(LEN=*) , INTENT(in ) :: clvar ! variable name 635 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid639 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional) 636 640 INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice) 637 641 INTEGER, DIMENSION(:) , INTENT(in ) :: map ! global-to-local mapping indices 638 642 !! 639 INTEGER :: ipi ! length of boundary data on local process 640 INTEGER :: ipj ! length of dummy dimension ( = 1 ) 641 INTEGER :: ipk ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 642 INTEGER :: ilendta ! length of data in file 643 INTEGER :: idvar ! variable ID 644 INTEGER :: ib, ik ! loop counters 645 INTEGER :: ierr 646 !! 647 CHARACTER(len=80) :: zfile 643 INTEGER :: ipi ! length of boundary data on local process 644 INTEGER :: ipj ! length of dummy dimension ( = 1 ) 645 INTEGER :: ipk ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 646 INTEGER :: ilendta ! length of data in file 647 INTEGER :: idvar ! variable ID 648 INTEGER :: ib, ik ! loop counters 649 INTEGER :: ierr 650 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read ! work space for global data 648 651 !!--------------------------------------------------------------------- 649 652 653 #if defined key_obc 654 dta_read => dta_global 655 #endif 656 650 657 ipi = SIZE( dta, 1 ) 651 658 ipj = 1 … … 655 662 ilendta = iom_file(num)%dimsz(1,idvar) 656 663 IF(lwp) WRITE(numout,*) 'Dim size for ',TRIM(clvar),' is ', ilendta 657 658 CALL iom_get ( num, jpdom_unknown, clvar, dta_global(1:ilendta,1:ipj,1:ipk), nrec ) 664 IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk 665 666 SELECT CASE( ipk ) 667 CASE(1) 668 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1 ), nrec ) 669 CASE DEFAULT 670 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec ) 671 END SELECT 659 672 ! 660 673 DO ib = 1, ipi 661 674 DO ik = 1, ipk 662 dta(ib,1,ik) = dta_ global(map(ib),1,ik)675 dta(ib,1,ik) = dta_read(map(ib),1,ik) 663 676 END DO 664 677 END DO 665 678 666 679 END SUBROUTINE fld_map 680 667 681 668 682 SUBROUTINE fld_rot( kt, sd ) -
branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r2797 r2814 47 47 USE obcini ! open boundary cond. initialization (obc_init routine) 48 48 USE obcdta ! open boundary cond. initialization (obc_dta_init routine) 49 USE obctides ! open boundary cond. initialization (tide_init routine) 49 50 USE istate ! initial state setting (istate_init routine) 50 51 USE ldfdyn ! lateral viscosity setting (ldfdyn_init routine) … … 296 297 IF( lk_obc ) CALL obc_init ! Open boundaries initialisation 297 298 IF( lk_obc ) CALL obc_dta_init ! Open boundaries initialisation of external data arrays 299 IF( lk_obc ) CALL tide_init ! Open boundaries initialisation of tidal harmonic forcing 298 300 299 301 CALL istate_init ! ocean initial state (Dynamics and tracers)
Note: See TracChangeset
for help on using the changeset viewer.