Changeset 8059
- Timestamp:
- 2017-05-23T10:32:39+02:00 (7 years ago)
- Location:
- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 5 added
- 32 edited
-
BDY/bdy_oce.F90 (modified) (1 diff)
-
BDY/bdydta.F90 (modified) (1 diff)
-
BDY/bdydyn3d.F90 (modified) (3 diffs)
-
BDY/bdyini.F90 (modified) (2 diffs)
-
BDY/bdylib.F90 (modified) (3 diffs)
-
BDY/bdytides.F90 (modified) (4 diffs)
-
BDY/bdytra.F90 (modified) (4 diffs)
-
DIA/dia25h.F90 (added)
-
DIA/diaharmana.F90 (added)
-
DIA/diainsitutem.F90 (added)
-
DIA/diatmb.F90 (added)
-
DIA/diawri.F90 (modified) (4 diffs)
-
DOM/domain.F90 (modified) (2 diffs)
-
DOM/dommsk.F90 (modified) (1 diff)
-
DOM/dtatsd.F90 (modified) (1 diff)
-
DYN/dynkeg.F90 (modified) (4 diffs)
-
DYN/dynspg_ts.F90 (modified) (4 diffs)
-
ICB/icbrst.F90 (modified) (3 diffs)
-
IOM/in_out_manager.F90 (modified) (1 diff)
-
IOM/restart.F90 (modified) (3 diffs)
-
SBC/sbc_oce.F90 (modified) (2 diffs)
-
SBC/sbcflx.F90 (modified) (7 diffs)
-
SBC/sbcssr.F90 (modified) (5 diffs)
-
SBC/tide.h90 (modified) (1 diff)
-
SBC/tide_mod.F90 (modified) (1 diff)
-
TRA/eosbn2.F90 (modified) (1 diff)
-
TRA/traadv_tvd.F90 (modified) (1 diff)
-
TRA/tradwl.F90 (added)
-
TRA/tranxt.F90 (modified) (3 diffs)
-
TRA/traqsr.F90 (modified) (9 diffs)
-
TRA/trasbc.F90 (modified) (2 diffs)
-
TRD/trdtra.F90 (modified) (1 diff)
-
ZDF/zdfmxl.F90 (modified) (7 diffs)
-
nemogcm.F90 (modified) (4 diffs)
-
step.F90 (modified) (4 diffs)
-
step_oce.F90 (modified) (3 diffs)
-
trc_oce.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
r8058 r8059 71 71 REAL, POINTER, DIMENSION(:,:) :: ht_s !: now snow thickness 72 72 #endif 73 #if defined key_top 74 CHARACTER(LEN=20) :: cn_obc !: type of boundary condition to apply 75 REAL(wp) :: rn_fac !: multiplicative scaling factor 76 REAL(wp), POINTER, DIMENSION(:,:) :: trc !: now field of the tracer 77 LOGICAL :: dmp !: obc damping term 78 #endif 79 73 80 END TYPE OBC_DATA 74 81 -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
r8058 r8059 37 37 #endif 38 38 USE sbcapr 39 #if defined key_top 40 USE par_trc 41 USE trc, ONLY: trn 42 #endif 39 43 40 44 IMPLICIT NONE -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90
r8058 r8059 61 61 CASE('zero') 62 62 CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 63 CASE('zerograd') 64 CALL bdy_dyn3d_zgrad( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 65 CASE('neumann') 66 CALL bdy_dyn3d_nmn( idx_bdy(ib_bdy), ib_bdy ) 63 67 CASE('orlanski') 64 68 CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. ) … … 117 121 118 122 END SUBROUTINE bdy_dyn3d_spe 123 124 SUBROUTINE bdy_dyn3d_zgrad( idx, dta, kt , ib_bdy ) 125 !!---------------------------------------------------------------------- 126 !! *** SUBROUTINE bdy_dyn3d_zgrad *** 127 !! 128 !! ** Purpose : - Enforce a zero gradient of normal velocity 129 !! 130 !!---------------------------------------------------------------------- 131 INTEGER :: kt 132 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 133 TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data 134 INTEGER, INTENT(in) :: ib_bdy ! BDY set index 135 !! 136 INTEGER :: jb, jk ! dummy loop indices 137 INTEGER :: ii, ij, igrd ! local integers 138 REAL(wp) :: zwgt ! boundary weight 139 INTEGER :: fu, fv 140 !!---------------------------------------------------------------------- 141 ! 142 IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zgrad') 143 ! 144 igrd = 2 ! Copying tangential velocity into bdy points 145 DO jb = 1, idx%nblenrim(igrd) 146 DO jk = 1, jpkm1 147 ii = idx%nbi(jb,igrd) 148 ij = idx%nbj(jb,igrd) 149 fu = ABS( ABS (NINT( idx%flagu(jb,igrd) ) ) - 1 ) 150 ua(ii,ij,jk) = ua(ii,ij,jk) * REAL( 1 - fu) + ( ua(ii,ij+fu,jk) * umask(ii,ij+fu,jk) & 151 &+ ua(ii,ij-fu,jk) * umask(ii,ij-fu,jk) ) * umask(ii,ij,jk) * REAL( fu ) 152 END DO 153 END DO 154 ! 155 igrd = 3 ! Copying tangential velocity into bdy points 156 DO jb = 1, idx%nblenrim(igrd) 157 DO jk = 1, jpkm1 158 ii = idx%nbi(jb,igrd) 159 ij = idx%nbj(jb,igrd) 160 fv = ABS( ABS (NINT( idx%flagv(jb,igrd) ) ) - 1 ) 161 va(ii,ij,jk) = va(ii,ij,jk) * REAL( 1 - fv ) + ( va(ii+fv,ij,jk) * vmask(ii+fv,ij,jk) & 162 &+ va(ii-fv,ij,jk) * vmask(ii-fv,ij,jk) ) * vmask(ii,ij,jk) * REAL( fv ) 163 END DO 164 END DO 165 CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy ) ! Boundary points should be updated 166 CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy ) 167 ! 168 IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 169 170 IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zgrad') 171 172 END SUBROUTINE bdy_dyn3d_zgrad 119 173 120 174 SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy ) … … 303 357 END SUBROUTINE bdy_dyn3d_dmp 304 358 359 SUBROUTINE bdy_dyn3d_nmn( idx, ib_bdy ) 360 !!---------------------------------------------------------------------- 361 !! *** SUBROUTINE bdy_dyn3d_nmn *** 362 !! 363 !! - Apply Neumann condition to baroclinic velocities. 364 !! - Wrapper routine for bdy_nmn 365 !! 366 !! 367 !!---------------------------------------------------------------------- 368 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 369 INTEGER, INTENT(in) :: ib_bdy ! BDY set index 370 371 INTEGER :: jb, igrd ! dummy loop indices 372 !!---------------------------------------------------------------------- 373 374 IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_nmn') 375 ! 376 !! Note that at this stage the ub and ua arrays contain the baroclinic velocities. 377 ! 378 igrd = 2 ! Neumann bc on u-velocity; 379 ! 380 CALL bdy_nmn( idx, igrd, ua ) 381 382 igrd = 3 ! Neumann bc on v-velocity 383 ! 384 CALL bdy_nmn( idx, igrd, va ) 385 ! 386 CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy ) ! Boundary points should be updated 387 CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy ) 388 ! 389 IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_nmn') 390 ! 391 END SUBROUTINE bdy_dyn3d_nmn 392 393 305 394 #else 306 395 !!---------------------------------------------------------------------- -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
r8058 r8059 213 213 dta_bdy(ib_bdy)%ll_u3d = .true. 214 214 dta_bdy(ib_bdy)%ll_v3d = .true. 215 CASE('neumann') 216 IF(lwp) WRITE(numout,*) ' Neumann conditions' 217 dta_bdy(ib_bdy)%ll_u3d = .false. 218 dta_bdy(ib_bdy)%ll_v3d = .false. 219 CASE('zerograd') 220 IF(lwp) WRITE(numout,*) ' Zero gradient for baroclinic velocities' 221 dta_bdy(ib_bdy)%ll_u3d = .false. 222 dta_bdy(ib_bdy)%ll_v3d = .false. 215 223 CASE('zero') 216 224 IF(lwp) WRITE(numout,*) ' Zero baroclinic velocities (runoff case)' … … 1087 1095 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 1088 1096 nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 1089 idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) ! tanh formulation 1097 idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) * 0.5 & 1098 & *(10./ FLOAT(nn_rimwidth(ib_bdy))) ) ! JGraham:modified for rim=15 1099 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) ! tanh formulation 1090 1100 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2. ! quadratic 1091 1101 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)) ! linear -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90
r8058 r8059 26 26 PUBLIC bdy_orlanski_2d ! routine called where? 27 27 PUBLIC bdy_orlanski_3d ! routine called where? 28 PUBLIC bdy_nmn ! routine called where? 28 29 29 30 !!---------------------------------------------------------------------- … … 354 355 END SUBROUTINE bdy_orlanski_3d 355 356 357 SUBROUTINE bdy_nmn( idx, igrd, phia ) 358 !!---------------------------------------------------------------------- 359 !! *** SUBROUTINE bdy_nmn *** 360 !! 361 !! ** Purpose : Duplicate the value at open boundaries, zero gradient. 362 !! 363 !!---------------------------------------------------------------------- 364 INTEGER, INTENT(in) :: igrd ! grid index 365 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: phia ! model after 3D field (to be updated) 366 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 367 !! 368 REAL(wp) :: zcoef, zcoef1, zcoef2 369 REAL(wp), POINTER, DIMENSION(:,:,:) :: pmask ! land/sea mask for field 370 REAL(wp), POINTER, DIMENSION(:,:) :: bdypmask ! land/sea mask for field 371 INTEGER :: ib, ik ! dummy loop indices 372 INTEGER :: ii, ij, ip, jp ! 2D addresses 373 !!---------------------------------------------------------------------- 374 ! 375 IF( nn_timing == 1 ) CALL timing_start('bdy_nmn') 376 ! 377 SELECT CASE(igrd) 378 CASE(1) 379 pmask => tmask(:,:,:) 380 bdypmask => bdytmask(:,:) 381 CASE(2) 382 pmask => umask(:,:,:) 383 bdypmask => bdyumask(:,:) 384 CASE(3) 385 pmask => vmask(:,:,:) 386 bdypmask => bdyvmask(:,:) 387 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for igrd in bdy_nmn' ) 388 END SELECT 389 DO ib = 1, idx%nblenrim(igrd) 390 ii = idx%nbi(ib,igrd) 391 ij = idx%nbj(ib,igrd) 392 DO ik = 1, jpkm1 393 ! search the sense of the gradient 394 zcoef1 = bdypmask(ii-1,ij )*pmask(ii-1,ij,ik) + bdypmask(ii+1,ij )*pmask(ii+1,ij,ik) 395 zcoef2 = bdypmask(ii ,ij-1)*pmask(ii,ij-1,ik) + bdypmask(ii ,ij+1)*pmask(ii,ij+1,ik) 396 IF ( nint(zcoef1+zcoef2) == 0) THEN 397 ! corner **** we probably only want to set the tangentail component for the dynamics here 398 zcoef = pmask(ii-1,ij,ik) + pmask(ii+1,ij,ik) + pmask(ii,ij-1,ik) + pmask(ii,ij+1,ik) 399 IF (zcoef > .5_wp) THEN ! Only set none isolated points. 400 phia(ii,ij,ik) = phia(ii-1,ij ,ik) * pmask(ii-1,ij ,ik) + & 401 & phia(ii+1,ij ,ik) * pmask(ii+1,ij ,ik) + & 402 & phia(ii ,ij-1,ik) * pmask(ii ,ij-1,ik) + & 403 & phia(ii ,ij+1,ik) * pmask(ii ,ij+1,ik) 404 phia(ii,ij,ik) = ( phia(ii,ij,ik) / zcoef ) * pmask(ii,ij,ik) 405 ELSE 406 phia(ii,ij,ik) = phia(ii,ij ,ik) * pmask(ii,ij ,ik) 407 ENDIF 408 ELSEIF ( nint(zcoef1+zcoef2) == 2) THEN 409 ! oblique corner **** we probably only want to set the normal component for the dynamics here 410 zcoef = pmask(ii-1,ij,ik)*bdypmask(ii-1,ij ) + pmask(ii+1,ij,ik)*bdypmask(ii+1,ij ) + & 411 & pmask(ii,ij-1,ik)*bdypmask(ii,ij -1 ) + pmask(ii,ij+1,ik)*bdypmask(ii,ij+1 ) 412 phia(ii,ij,ik) = phia(ii-1,ij ,ik) * pmask(ii-1,ij ,ik)*bdypmask(ii-1,ij ) + & 413 & phia(ii+1,ij ,ik) * pmask(ii+1,ij ,ik)*bdypmask(ii+1,ij ) + & 414 & phia(ii ,ij-1,ik) * pmask(ii ,ij-1,ik)*bdypmask(ii,ij -1 ) + & 415 & phia(ii ,ij+1,ik) * pmask(ii ,ij+1,ik)*bdypmask(ii,ij+1 ) 416 417 phia(ii,ij,ik) = ( phia(ii,ij,ik) / MAX(1._wp, zcoef) ) * pmask(ii,ij,ik) 418 ELSE 419 ip = nint(bdypmask(ii+1,ij )*pmask(ii+1,ij,ik) - bdypmask(ii-1,ij )*pmask(ii-1,ij,ik)) 420 jp = nint(bdypmask(ii ,ij+1)*pmask(ii,ij+1,ik) - bdypmask(ii ,ij-1)*pmask(ii,ij-1,ik)) 421 phia(ii,ij,ik) = phia(ii+ip,ij+jp,ik) * pmask(ii+ip,ij+jp,ik) 422 ENDIF 423 END DO 424 END DO 425 ! 426 IF( nn_timing == 1 ) CALL timing_stop('bdy_nmn') 427 ! 428 END SUBROUTINE bdy_nmn 429 356 430 357 431 #else … … 366 440 WRITE(*,*) 'bdy_orlanski_3d: You should not have seen this print! error?', kt 367 441 END SUBROUTINE bdy_orlanski_3d 442 SUBROUTINE bdy_nmn( idx, igrd, phia ) ! Empty routine 443 WRITE(*,*) 'bdy_nmn: You should not have seen this print! error?', kt 444 END SUBROUTINE bdy_nmn 368 445 #endif 369 446 -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
r8058 r8059 50 50 REAL(wp), POINTER, DIMENSION(:,:,:) :: v !: Tidal constituents : V (after nodal cor.) 51 51 END TYPE TIDES_DATA 52 INTEGER, PUBLIC, PARAMETER :: jptides_max = 15 !: Max number of tidal contituents 53 LOGICAL, PUBLIC :: ln_harm_ana_store !: =T Stores data for harmonic Analysis 54 LOGICAL, PUBLIC :: ln_harm_ana_compute !: =T Compute harmonic Analysis 55 LOGICAL, PUBLIC :: ln_harmana_read !: =T Decide to do the analysis 56 !from scratch or continue previous run 52 57 53 58 !$AGRIF_DO_NOT_TREAT … … 90 95 TYPE(MAP_POINTER), DIMENSION(jpbgrd) :: ibmap_ptr !: array of pointers to nbmap 91 96 !! 92 NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj 97 NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj, ln_harm_ana_store, ln_harm_ana_compute, ln_harmana_read 93 98 !!---------------------------------------------------------------------- 94 99 … … 102 107 103 108 REWIND(numnam_cfg) 109 REWIND(numnam_ref) ! slwa 104 110 105 111 DO ib_bdy = 1, nb_bdy … … 125 131 IF(lwp) WRITE(numout,*) ' assume complex conjugate : ', ln_bdytide_conj 126 132 IF(lwp) WRITE(numout,*) ' Number of tidal components to read: ', nb_harmo 133 IF(lwp) WRITE(numout,*) ' Use PCOMS harmonic ananalysis or not: ', ln_harm_ana_store 134 IF(lwp) WRITE(numout,*) ' Compute Final harmonic ananalysis or not: ', ln_harm_ana_compute 135 IF(lwp) WRITE(numout,*) ' Read in previous days harmonic data or start afresh: ', ln_harmana_read 127 136 IF(lwp) THEN 128 137 WRITE(numout,*) ' Tidal components: ' -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90
r8058 r8059 91 91 !! 92 92 REAL(wp) :: zwgt ! boundary weight 93 INTEGER :: ib, ik, igrd ! dummy loop indices 94 INTEGER :: ii, ij ! 2D addresses 93 REAL(wp) :: zcoef, zcoef1,zcoef2 94 INTEGER :: ib, ik, igrd ! dummy loop indices 95 INTEGER :: ii, ij, ip, jp ! 2D addresses 95 96 !!---------------------------------------------------------------------- 96 97 ! … … 160 161 !! 161 162 REAL(wp) :: zwgt ! boundary weight 162 INTEGER :: ib, ik, igrd ! dummy loop indices 163 INTEGER :: ii, ij,zcoef, zcoef1,zcoef2, ip, jp ! 2D addresses 163 REAL(wp) :: zcoef, zcoef1,zcoef2 164 INTEGER :: ib, ik, igrd ! dummy loop indices 165 INTEGER :: ii, ij, ip, jp ! 2D addresses 164 166 !!---------------------------------------------------------------------- 165 167 ! … … 174 176 zcoef1 = bdytmask(ii-1,ij ) + bdytmask(ii+1,ij ) 175 177 zcoef2 = bdytmask(ii ,ij-1) + bdytmask(ii ,ij+1) 176 IF ( zcoef1+zcoef2 == 0) THEN178 IF ( NINT(zcoef1+zcoef2) == 0) THEN 177 179 ! corner 178 180 zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) + tmask(ii,ij-1,ik) + tmask(ii,ij+1,ik) … … 181 183 & tsa(ii ,ij-1,ik,jp_tem) * tmask(ii ,ij-1,ik) + & 182 184 & tsa(ii ,ij+1,ik,jp_tem) * tmask(ii ,ij+1,ik) 183 tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1 , zcoef) ) * tmask(ii,ij,ik)185 tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1._wp, zcoef) ) * tmask(ii,ij,ik) 184 186 tsa(ii,ij,ik,jp_sal) = tsa(ii-1,ij ,ik,jp_sal) * tmask(ii-1,ij ,ik) + & 185 187 & tsa(ii+1,ij ,ik,jp_sal) * tmask(ii+1,ij ,ik) + & 186 188 & tsa(ii ,ij-1,ik,jp_sal) * tmask(ii ,ij-1,ik) + & 187 189 & tsa(ii ,ij+1,ik,jp_sal) * tmask(ii ,ij+1,ik) 188 tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1 , zcoef) ) * tmask(ii,ij,ik)190 tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1._wp, zcoef) ) * tmask(ii,ij,ik) 189 191 ELSE 190 ip = bdytmask(ii+1,ij ) - bdytmask(ii-1,ij )191 jp = bdytmask(ii ,ij+1) - bdytmask(ii ,ij-1)192 ip = NINT(bdytmask(ii+1,ij ) - bdytmask(ii-1,ij )) 193 jp = NINT(bdytmask(ii ,ij+1) - bdytmask(ii ,ij-1)) 192 194 tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii+ip,ij+jp,ik) 193 195 tsa(ii,ij,ik,jp_sal) = tsa(ii+ip,ij+jp,ik,jp_sal) * tmask(ii+ip,ij+jp,ik) -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
r8058 r8059 44 44 USE in_out_manager ! I/O manager 45 45 USE diadimg ! dimg direct access file format output 46 USE diatmb ! Top,middle,bottom output 47 USE dia25h ! 25h Mean output 46 48 USE iom 47 49 USE ioipsl 48 50 USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities 51 USE eosbn2 ! equation of state (eos_bn2 routine) 52 49 53 50 54 #if defined key_lim2 … … 132 136 REAL(wp), POINTER, DIMENSION(:,:) :: z2d ! 2D workspace 133 137 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d ! 3D workspace 138 REAL(wp), POINTER, DIMENSION(:,:,:) :: zrhd , zrhop ! 3D workspace 134 139 !!---------------------------------------------------------------------- 135 140 ! … … 138 143 CALL wrk_alloc( jpi , jpj , z2d ) 139 144 CALL wrk_alloc( jpi , jpj, jpk , z3d ) 145 CALL wrk_alloc( jpi , jpj, jpk , zrhd , zrhop ) 140 146 ! 141 147 ! Output the initial state and forcings … … 376 382 CALL iom_put( "v_salttr", 0.5 * z2d ) ! heat transport in j-direction 377 383 ENDIF 384 385 IF( iom_use("rhop") ) THEN 386 CALL eos( tsn, zrhd, zrhop, fsdept_n(:,:,:) ) ! now in situ and potential density 387 zrhop(:,:,jpk) = 0._wp 388 CALL iom_put( 'rhop', zrhop ) 389 ENDIF 390 378 391 ! 379 392 CALL wrk_dealloc( jpi , jpj , z2d ) 380 393 CALL wrk_dealloc( jpi , jpj, jpk , z3d ) 394 CALL wrk_dealloc( jpi , jpj, jpk , zrhd , zrhop ) 395 ! 396 ! If we want tmb values 397 398 IF (ln_diatmb) THEN 399 CALL dia_tmb 400 ENDIF 401 IF (ln_dia25h) THEN 402 CALL dia_25h( kt ) 403 ENDIF 381 404 ! 382 405 IF( nn_timing == 1 ) CALL timing_stop('dia_wri') -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r8058 r8059 136 136 USE ioipsl 137 137 NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list, & 138 & nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rst art , nn_rstctl, &138 & nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstdate, ln_rstart , nn_rstctl, & 139 139 & nn_it000, nn_itend , nn_date0 , nn_leapy , nn_istate , nn_stock , & 140 140 & nn_write, ln_dimgnnn, ln_mskland , ln_cfmeta , ln_clobber, nn_chunksz, nn_euler … … 174 174 WRITE(numout,*) ' restart output directory cn_ocerst_outdir= ', cn_ocerst_outdir 175 175 WRITE(numout,*) ' restart logical ln_rstart = ', ln_rstart 176 WRITE(numout,*) ' datestamping of restarts ln_rstdate = ', ln_rstdate 176 177 WRITE(numout,*) ' start with forward time step nn_euler = ', nn_euler 177 178 WRITE(numout,*) ' control of time step nn_rstctl = ', nn_rstctl -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r8058 r8059 31 31 USE wrk_nemo ! Memory allocation 32 32 USE timing ! Timing 33 USE iom ! slwa 33 34 34 35 IMPLICIT NONE -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
r8058 r8059 209 209 ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:) 210 210 ! 211 IF( ln_sco ) THEN !== s- or mixed s-zps-coordinate ==!211 IF( ln_sco .and. 1==0) THEN !== s- or mixed s-zps-coordinate ==! 212 212 ! 213 213 CALL wrk_alloc( jpk, ztp, zsp ) -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
r8058 r8059 24 24 USE wrk_nemo ! Memory Allocation 25 25 USE timing ! Timing 26 #if defined key_bdy 27 USE bdy_oce ! ocean open boundary conditions 28 #endif 26 29 27 30 IMPLICIT NONE … … 78 81 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke 79 82 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 83 #if defined key_bdy 84 INTEGER :: jb ! dummy loop indices 85 INTEGER :: ii, ij, igrd, ib_bdy ! local integers 86 INTEGER :: fu, fv 87 #endif 80 88 !!---------------------------------------------------------------------- 81 89 ! … … 97 105 98 106 zhke(:,:,jpk) = 0._wp 99 107 108 #if defined key_bdy 109 ! Maria Luneva & Fred Wobus: July-2016 110 ! compensate for lack of turbulent kinetic energy on liquid bdy points 111 DO ib_bdy = 1, nb_bdy 112 IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 113 igrd = 2 ! Copying normal velocity into points outside bdy 114 DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 115 DO jk = 1, jpkm1 116 ii = idx_bdy(ib_bdy)%nbi(jb,igrd) 117 ij = idx_bdy(ib_bdy)%nbj(jb,igrd) 118 fu = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) ) 119 un(ii-fu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk) 120 END DO 121 END DO 122 ! 123 igrd = 3 ! Copying normal velocity into points outside bdy 124 DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 125 DO jk = 1, jpkm1 126 ii = idx_bdy(ib_bdy)%nbi(jb,igrd) 127 ij = idx_bdy(ib_bdy)%nbj(jb,igrd) 128 fv = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) ) 129 vn(ii,ij-fv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk) 130 END DO 131 END DO 132 ENDIF 133 ENDDO 134 #endif 135 100 136 SELECT CASE ( kscheme ) !== Horizontal kinetic energy at T-point ==! 101 137 ! … … 134 170 END SELECT 135 171 ! 172 173 #if defined key_bdy 174 ! restore velocity masks at points outside boundary 175 un(:,:,:) = un(:,:,:) * umask(:,:,:) 176 vn(:,:,:) = vn(:,:,:) * vmask(:,:,:) 177 #endif 178 136 179 DO jk = 1, jpkm1 !== grad( KE ) added to the general momentum trends ==! 137 180 DO jj = 2, jpjm1 -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r8058 r8059 41 41 USE timing ! Timing 42 42 USE sbcapr ! surface boundary condition: atmospheric pressure 43 USE diatmb ! Top,middle,bottom output 43 44 USE dynadv, ONLY: ln_dynadv_vec 44 45 #if defined key_agrif … … 144 145 INTEGER :: ji, jj, jk, jn ! dummy loop indices 145 146 INTEGER :: ikbu, ikbv, noffset ! local integers 147 REAL(wp) :: zmdi 146 148 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 147 149 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - … … 169 171 CALL wrk_alloc( jpi, jpj, zhf ) 170 172 ! 173 zmdi=1.e+20 ! missing data indicator for masking 171 174 ! !* Local constant initialization 172 175 z1_12 = 1._wp / 12._wp … … 926 929 CALL wrk_dealloc( jpi, jpj, zhf ) 927 930 ! 931 IF ( ln_diatmb ) THEN 932 CALL iom_put( "baro_u" , un_b*umask(:,:,1)+zmdi*(1-umask(:,:,1 ) ) ) ! Barotropic U Velocity 933 CALL iom_put( "baro_v" , vn_b*vmask(:,:,1)+zmdi*(1-vmask(:,:,1 ) ) ) ! Barotropic V Velocity 934 ENDIF 928 935 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') 929 936 ! -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90
r8058 r8059 18 18 !!---------------------------------------------------------------------- 19 19 USE par_oce ! NEMO parameters 20 USE phycst ! for rday 20 21 USE dom_oce ! NEMO domain 21 22 USE in_out_manager ! NEMO IO routines 23 USE ioipsl, ONLY : ju2ymds ! for calendar 22 24 USE lib_mpp ! NEMO MPI library, lk_mpp in particular 23 25 USE netcdf ! netcdf routines for IO … … 231 233 INTEGER :: jn ! dummy loop index 232 234 INTEGER :: ix_dim, iy_dim, ik_dim, in_dim 233 CHARACTER(len=256) :: cl_path 234 CHARACTER(len=256) :: cl_filename 235 INTEGER :: iyear, imonth, iday 236 REAL (wp) :: zsec 237 REAL (wp) :: zfjulday 238 CHARACTER(len=256) :: cl_path 239 CHARACTER(len=256) :: cl_filename 240 CHARACTER(LEN=20) :: clkt ! ocean time-step deine as a character 235 241 TYPE(iceberg), POINTER :: this 236 242 TYPE(point) , POINTER :: pt … … 240 246 cl_path = TRIM(cn_ocerst_outdir) 241 247 IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/' 248 IF ( ln_rstdate ) THEN 249 zfjulday = fjulday + rdttra(1) / rday 250 IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday ) zfjulday = REAL(NINT(zfjulday),wp) ! avoid truncation error 251 CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec ) 252 WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 253 ELSE 254 IF( kt > 999999999 ) THEN ; WRITE(clkt, * ) kt 255 ELSE ; WRITE(clkt, '(i8.8)') kt 256 ENDIF 257 ENDIF 242 258 IF( lk_mpp ) THEN 243 WRITE(cl_filename,'(A,"_icebergs_", I8.8,"_restart_",I4.4,".nc")') TRIM(cexper), kt, narea-1259 WRITE(cl_filename,'(A,"_icebergs_",A,"_restart_",I4.4,".nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)), narea-1 244 260 ELSE 245 WRITE(cl_filename,'(A,"_icebergs_", I8.8,"_restart.nc")') TRIM(cexper), kt261 WRITE(cl_filename,'(A,"_icebergs_",A,"_restart.nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)) 246 262 ENDIF 247 263 IF (nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, write_restart: creating ',TRIM(cl_path)//TRIM(cl_filename) -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90
r8058 r8059 30 30 CHARACTER(lc) :: cn_ocerst_outdir !: restart output directory 31 31 LOGICAL :: ln_rstart !: start from (F) rest or (T) a restart file 32 LOGICAL :: ln_rstdate !: datestamping of restarts 32 33 LOGICAL :: ln_rst_list !: output restarts at list of times (T) or by frequency (F) 33 34 INTEGER :: nn_no !: job number -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
r8058 r8059 21 21 USE in_out_manager ! I/O manager 22 22 USE iom ! I/O module 23 USE ioipsl, ONLY : ju2ymds ! for calendar 23 24 USE eosbn2 ! equation of state (eos bn2 routine) 24 25 USE trdmxl_oce ! ocean active mixed layer tracers trends variables … … 54 55 !!---------------------------------------------------------------------- 55 56 INTEGER, INTENT(in) :: kt ! ocean time-step 57 INTEGER :: iyear, imonth, iday 58 REAL (wp) :: zsec 59 REAL (wp) :: zfjulday 56 60 !! 57 61 CHARACTER(LEN=20) :: clkt ! ocean time-step deine as a character 58 62 CHARACTER(LEN=50) :: clname ! ocean output restart file name 59 CHARACTER( lc):: clpath ! full path to ocean output restart file63 CHARACTER(LEN=150) :: clpath ! full path to ocean output restart file 60 64 !!---------------------------------------------------------------------- 61 65 ! … … 81 85 IF( kt == nitrst - 1 .OR. nstock == 1 .OR. ( kt == nitend .AND. .NOT. lrst_oce ) ) THEN 82 86 IF( nitrst <= nitend .AND. nitrst > 0 ) THEN 83 ! beware of the format used to write kt (default is i8.8, that should be large enough...) 84 IF( nitrst > 999999999 ) THEN ; WRITE(clkt, * ) nitrst 85 ELSE ; WRITE(clkt, '(i8.8)') nitrst 87 IF ( ln_rstdate ) THEN 88 zfjulday = fjulday + rdttra(1) / rday 89 IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday ) zfjulday = REAL(NINT(zfjulday),wp) ! avoid truncation error 90 CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec ) 91 WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 92 ELSE 93 ! beware of the format used to write kt (default is i8.8, that should be large enough...) 94 IF( nitrst > 999999999 ) THEN ; WRITE(clkt, * ) nitrst 95 ELSE ; WRITE(clkt, '(i8.8)') nitrst 96 ENDIF 86 97 ENDIF 87 98 ! create the file -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
r8058 r8059 121 121 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sprecip !: solid precipitation [Kg/m2/s] 122 122 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fr_i !: ice fraction = 1 - lead fraction (between 0 to 1) 123 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: pressnow !: UKMO SHELF pressure 124 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: apgu !: UKMO SHELF pressure forcing 125 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: apgv !: UKMO SHELF pressure forcing 123 126 #if defined key_cpl_carbon_cycle 124 127 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: atm_co2 !: atmospheric pCO2 [ppm] … … 171 174 #endif 172 175 & ssu_m (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) , & 176 & pressnow(jpi,jpj), apgu(jpi,jpj) , apgv(jpi,jpj) , & 173 177 & ssv_m (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) ) 174 178 ! -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90
r8058 r8059 28 28 PUBLIC sbc_flx ! routine called by step.F90 29 29 30 INTEGER , PARAMETER :: jpfld = 5! maximum number of files to read30 INTEGER , PARAMETER :: jpfld = 6 ! maximum number of files to read 31 31 INTEGER , PARAMETER :: jp_utau = 1 ! index of wind stress (i-component) file 32 32 INTEGER , PARAMETER :: jp_vtau = 2 ! index of wind stress (j-component) file … … 34 34 INTEGER , PARAMETER :: jp_qsr = 4 ! index of solar heat file 35 35 INTEGER , PARAMETER :: jp_emp = 5 ! index of evaporation-precipation file 36 INTEGER , PARAMETER :: jp_press = 6 ! index of pressure for UKMO shelf fluxes 36 37 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read) 38 LOGICAL , PUBLIC :: ln_shelf_flx = .FALSE. ! UKMO SHELF specific flux flag 39 INTEGER :: jpfld_local ! maximum number of files to read (locally modified depending on ln_shelf_flx) 37 40 38 41 !! * Substitutions … … 82 85 REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient 83 86 REAL(wp) :: ztx, zty, zmod, zcoef ! temporary variables 87 REAL :: cs ! UKMO SHELF: Friction co-efficient at surface 88 REAL :: totwindspd ! UKMO SHELF: Magnitude of wind speed vector 89 90 REAL(wp) :: rhoa = 1.22 ! Air density kg/m3 91 REAL(wp) :: cdrag = 1.5e-3 ! drag coefficient 84 92 !! 85 93 CHARACTER(len=100) :: cn_dir ! Root directory for location of flx files 86 94 TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist information structures 87 TYPE(FLD_N) :: sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp ! informations about the fields to be read 88 NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp 95 TYPE(FLD_N) :: sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp, sn_press ! informations about the fields to be read 96 LOGICAL :: ln_foam_flx = .FALSE. ! UKMO FOAM specific flux flag 97 NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp, & 98 & ln_foam_flx, sn_press, ln_shelf_flx 89 99 !!--------------------------------------------------------------------- 90 100 ! … … 109 119 slf_i(jp_emp ) = sn_emp 110 120 ! 111 ALLOCATE( sf(jpfld), STAT=ierror ) ! set sf structure 121 ALLOCATE( sf(jpfld), STAT=ierror ) ! set sf structure 122 IF( ln_shelf_flx ) slf_i(jp_press) = sn_press 123 124 ! define local jpfld depending on shelf_flx logical 125 IF( ln_shelf_flx ) THEN 126 jpfld_local = jpfld 127 ELSE 128 jpfld_local = jpfld-1 129 ENDIF 130 ! 112 131 IF( ierror > 0 ) THEN 113 132 CALL ctl_stop( 'sbc_flx: unable to allocate sf structure' ) ; RETURN 114 133 ENDIF 115 DO ji= 1, jpfld 134 DO ji= 1, jpfld_local 116 135 ALLOCATE( sf(ji)%fnow(jpi,jpj,1) ) 117 136 IF( slf_i(ji)%ln_tint ) ALLOCATE( sf(ji)%fdta(jpi,jpj,1,2) ) … … 132 151 ENDIF 133 152 !CDIR COLLAPSE 153 !!UKMO SHELF effect of atmospheric pressure on SSH 154 ! If using ln_apr_dyn, this is done there so don't repeat here. 155 IF( ln_shelf_flx .AND. .NOT. ln_apr_dyn) THEN 156 DO jj = 1, jpjm1 157 DO ji = 1, jpim1 158 apgu(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji+1,jj,1)-sf(jp_press)%fnow(ji,jj,1))/e1u(ji,jj) 159 apgv(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji,jj+1,1)-sf(jp_press)%fnow(ji,jj,1))/e2v(ji,jj) 160 END DO 161 END DO 162 ENDIF ! ln_shelf_flx 163 134 164 DO jj = 1, jpj ! set the ocean fluxes from read fields 135 165 DO ji = 1, jpi 136 utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) 137 vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) 138 qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) 139 emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1) 166 IF( ln_shelf_flx ) THEN 167 !! UKMO SHELF - need atmospheric pressure to calculate Haney forcing 168 pressnow(ji,jj) = sf(jp_press)%fnow(ji,jj,1) 169 !! UKMO SHELF flux files contain wind speed not wind stress 170 totwindspd = sqrt((sf(jp_utau)%fnow(ji,jj,1))**2.0 + (sf(jp_vtau)%fnow(ji,jj,1))**2.0) 171 cs = 0.63 + (0.066 * totwindspd) 172 utau(ji,jj) = cs * (rhoa/rau0) * sf(jp_utau)%fnow(ji,jj,1) * totwindspd 173 vtau(ji,jj) = cs * (rhoa/rau0) * sf(jp_vtau)%fnow(ji,jj,1) * totwindspd 174 ELSE 175 utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) 176 vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) 177 ENDIF 178 qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1) 179 IF( ln_foam_flx .OR. ln_shelf_flx ) THEN 180 !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot) 181 qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) 182 !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P 183 emp (ji,jj) = -1. * sf(jp_emp )%fnow(ji,jj,1) 184 ELSE 185 qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) 186 emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1) 187 ENDIF 140 188 END DO 141 189 END DO … … 143 191 qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp ! mass flux is at SST 144 192 ! 193 194 !! UKMO FOAM wind fluxes need lbc_lnk calls owing to a bug in interp.exe 195 IF( ln_foam_flx ) THEN 196 CALL lbc_lnk( utau(:,:), 'U', -1. ) 197 CALL lbc_lnk( vtau(:,:), 'V', -1. ) 198 ENDIF 199 145 200 ! ! module of wind stress and wind speed at T-point 146 201 zcoef = 1. / ( zrhoa * zcdrag ) … … 162 217 WRITE(numout,*) 163 218 WRITE(numout,*) ' read daily momentum, heat and freshwater fluxes OK' 164 DO jf = 1, jpfld 219 DO jf = 1, jpfld_local 165 220 IF( jf == jp_utau .OR. jf == jp_vtau ) zfact = 1. 166 221 IF( jf == jp_qtot .OR. jf == jp_qsr ) zfact = 0.1 -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90
r8058 r8059 42 42 LOGICAL :: ln_sssr_bnd ! flag to bound erp term 43 43 REAL(wp) :: rn_sssr_bnd ! ABS(Max./Min.) value of erp term [mm/day] 44 LOGICAL :: ln_UKMO_haney ! UKMO specific flag to calculate Haney forcing 44 45 45 46 REAL(wp) , ALLOCATABLE, DIMENSION(:) :: buffer ! Temporary buffer for exchange … … 79 80 INTEGER :: ierror ! return error code 80 81 !! 82 REAL(wp) :: sst1,sst2 ! sea surface temperature 83 REAL(wp) :: e_sst1, e_sst2 ! saturation vapour pressure 84 REAL(wp) :: qs1,qs2 ! specific humidity 85 REAL(wp) :: pr_tmp ! temporary variable for pressure 86 87 REAL(wp), DIMENSION(jpi,jpj) :: hny_frc1 ! Haney forcing for sensible heat, correction for latent heat 88 REAL(wp), DIMENSION(jpi,jpj) :: hny_frc2 ! Haney forcing for sensible heat, correction for latent heat 89 !! 81 90 CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files 82 91 TYPE(FLD_N) :: sn_sst, sn_sss ! informations about the fields to be read … … 95 104 ! 96 105 IF( nn_sstr == 1 ) THEN !* Temperature restoring term 97 DO jj = 1, jpj 98 DO ji = 1, jpi 99 zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) ) 100 qns(ji,jj) = qns(ji,jj) + zqrp 101 qrp(ji,jj) = zqrp 102 END DO 103 END DO 106 IF( ln_UKMO_haney ) THEN 107 DO jj = 1, jpj 108 DO ji = 1, jpi 109 sst1 = sst_m(ji,jj) 110 sst2 = sf_sst(1)%fnow(ji,jj,1) 111 e_sst1 = 10**((0.7859+0.03477*sst1)/(1.+0.00412*sst1)) 112 e_sst2 = 10**((0.7859+0.03477*sst2)/(1.+0.00412*sst2)) 113 pr_tmp = 0.01*pressnow(ji,jj) !pr_tmp = 1012.0 114 qs1 = (0.62197*e_sst1)/(pr_tmp-0.378*e_sst1) 115 qs2 = (0.62197*e_sst2)/(pr_tmp-0.378*e_sst2) 116 hny_frc1(ji,jj) = sst1-sst2 117 hny_frc2(ji,jj) = qs1-qs2 118 !Might need to mask off land points. 119 hny_frc1(ji,jj)=-hny_frc1(ji,jj)*wndm(ji,jj)*1.42 120 hny_frc2(ji,jj)=-hny_frc2(ji,jj)*wndm(ji,jj)*4688.0 121 qns(ji,jj)=qns(ji,jj)+hny_frc1(ji,jj)+hny_frc2(ji,jj) 122 qrp(ji,jj) = 0.e0 123 END DO 124 END DO 125 ELSE 126 DO jj = 1, jpj 127 DO ji = 1, jpi 128 zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) ) 129 qns(ji,jj) = qns(ji,jj) + zqrp 130 qrp(ji,jj) = zqrp 131 END DO 132 END DO 133 ENDIF 104 134 CALL iom_put( "qrp", qrp ) ! heat flux damping 105 135 ENDIF … … 163 193 CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files 164 194 TYPE(FLD_N) :: sn_sst, sn_sss ! informations about the fields to be read 165 NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, sn_sss, ln_sssr_bnd, rn_sssr_bnd 195 NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, sn_sss, ln_sssr_bnd, rn_sssr_bnd, ln_UKMO_haney 166 196 INTEGER :: ios 167 197 !!---------------------------------------------------------------------- … … 189 219 WRITE(numout,*) ' flag to bound erp term ln_sssr_bnd = ', ln_sssr_bnd 190 220 WRITE(numout,*) ' ABS(Max./Min.) erp threshold rn_sssr_bnd = ', rn_sssr_bnd, ' mm/day' 221 WRITE(numout,*) ' Haney forcing ln_UKMO_haney = ', ln_UKMO_haney 191 222 ENDIF 192 223 ! -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/tide.h90
r4292 r8059 28 28 Wave(18) = tide( 'L2' , 0.006694 , 2 , 2 , -1 , 2 , -1 , 0 , +180 , 2 , -2 , 0 , 0 , 0 , 215 ) 29 29 Wave(19) = tide( 'T2' , 0.006614 , 2 , 2 , 0 , -1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ) 30 ! 31 ! !! name_tide , equitide , nutide , nt , ns , nh , np , np1 , shift , nksi , nnu0 , nnu1 , nnu2 , R , formula !! 32 Wave(20) = tide( 'MNS2' , 0.000000 , 2 , 2 , -5 , 4 , 1 , 0 , 0 , 4 , -4 , 0 , 0 , 0 , 6 ) 33 Wave(21) = tide( 'Lam2' , 0.001760 , 2 , 2 , -1 , 0 , 1 , 0 , +180 , 2 , -2 , 0 , 0 , 0 , 78 ) 34 Wave(22) = tide( 'MSN2' , 0.000000 , 2 , 2 , 1 , 0 , 1 , 0 , 0 , 2 , -2 , 0 , 2 , 0 , 6 ) 35 Wave(23) = tide( '2SM2' , 0.000000 , 2 , 2 , 2 , -2 , 0 , 0 , 0 , -2 , 2 , 0 , 0 , 0 , 16 ) 36 Wave(24) = tide( 'MO3' , 0.000000 , 3 , 3 , -4 , 1 , 0 , 0 , +90 , 2 , -2 , 0 , 0 , 0 , 13 ) 37 Wave(25) = tide( 'MK3' , 0.000000 , 3 , 3 , -2 , 3 , 0 , 0 , -90 , 2 , -2 , -1 , 0 , 0 , 10 ) 38 Wave(26) = tide( 'MN4' , 0.000000 , 4 , 4 , -5 , 4 , 1 , 0 , 0 , 4 , -4 , 0 , 0 , 0 , 1 ) 39 Wave(27) = tide( 'MS4' , 0.000000 , 4 , 4 , -2 , 2 , 0 , 0 , 0 , 2 , -2 , 0 , 0 , 0 , 2 ) 40 Wave(28) = tide( 'M6' , 0.000000 , 6 , 6 , -6 , 6 , 0 , 0 , 0 , 6 , -6 , 0 , 0 , 0 , 4 ) 41 Wave(29) = tide( '2MS6' , 0.000000 , 6 , 6 , -4 , 4 , 0 , 0 , 0 , 4 , -4 , 0 , 0 , 0 , 6 ) 42 Wave(30) = tide( '2MK6' , 0.000000 , 6 , 6 , -4 , 6 , 0 , 0 , 0 , 4 , -4 , 0 , -2 , 0 , 5 ) 43 Wave(31) = tide( '3M2S2' , 0.000000 , 2 , 2 , -6 , 6 , 0 , 0 , 0 , 6 , -6 , 0 , 0 , 0 , 12 ) -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90
r8058 r8059 16 16 PUBLIC tide_init_Wave ! called by tideini and diaharm modules 17 17 18 INTEGER, PUBLIC, PARAMETER :: jpmax_harmo = 19!: maximum number of harmonic18 INTEGER, PUBLIC, PARAMETER :: jpmax_harmo = 31 !: maximum number of harmonic 19 19 20 20 TYPE, PUBLIC :: tide 21 CHARACTER(LEN= 4) :: cname_tide21 CHARACTER(LEN=5) :: cname_tide 22 22 REAL(wp) :: equitide 23 23 INTEGER :: nutide -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r8058 r8059 1241 1241 IF(lwm) WRITE( numond, nameos ) 1242 1242 ! 1243 rau0 = 1026._wp !: volumic mass of reference [kg/m3] 1243 rau0 = 1020._wp !: volumic mass of reference [kg/m3] 1244 ! rau0 = 1026._wp !: volumic mass of reference [kg/m3] 1244 1245 rcp = 3991.86795711963_wp !: heat capacity [J/K] 1245 1246 ! -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
r8058 r8059 100 100 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 101 101 ENDIF 102 ! slwa unless you use l_trdtra too, the above switches off trend calculations for l_trdtrc 103 l_trd = .FALSE. 104 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 105 !slwa 102 106 ! 103 107 IF( l_trd ) THEN -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
r8058 r8059 58 58 59 59 REAL(wp) :: rbcp ! Brown & Campana parameters for semi-implicit hpg 60 INTEGER :: warn_1, warn_2 ! indicators for warning statement 60 61 61 62 !! * Substitutions … … 93 94 INTEGER, INTENT(in) :: kt ! ocean time-step index 94 95 !! 95 INTEGER :: jk, jn ! dummy loop indices96 REAL(wp) :: zfact ! local scalars96 INTEGER :: jk, jn, ji, jj ! dummy loop indices 97 REAL(wp) :: zfact, zfreeze ! local scalars 97 98 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds 98 99 !!---------------------------------------------------------------------- … … 120 121 IF( lk_bdy ) CALL bdy_tra( kt ) ! BDY open boundaries 121 122 #endif 122 123 124 #if ( ! defined key_lim3 && ! defined key_lim2 && ! key_cice ) 125 IF ( kt == nit000 ) warn_1=0 126 warn_2=0 127 DO jk = 1, jpkm1 128 DO jj = 1, jpj 129 DO ji = 1, jpi 130 IF ( tsa(ji,jj,jk,jp_tem) .lt. 0.0 ) THEN 131 ! calculate freezing point 132 zfreeze = ( -0.0575_wp + 1.710523E-3 * Sqrt(Abs(tsn(ji,jj,jk,jp_sal))) & 133 - 2.154996E-4 * tsn(ji,jj,jk,jp_sal) ) * tsn(ji,jj,jk,jp_sal) - 7.53E-4 * ( 10.0_wp + fsdept(ji,jj,jk) ) 134 IF ( tsa(ji,jj,jk,jp_tem) .lt. zfreeze ) THEN 135 tsa(ji,jj,jk,jp_tem)=zfreeze 136 warn_2=1 137 ENDIF 138 ENDIF 139 END DO 140 END DO 141 END DO 142 CALL mpp_max(warn_1) 143 CALL mpp_max(warn_2) 144 IF ( (warn_1 == 0) .and. (warn_2 /= 0) ) THEN 145 IF(lwp) THEN 146 CALL ctl_warn( ' Temperatures dropping below freezing point, ', & 147 & ' being forced to freezing point, no longer conservative' ) 148 ENDIF 149 warn_1=1 150 ENDIF 151 #endif 152 123 153 ! set time step size (Euler/Leapfrog) 124 154 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dtra(:) = rdttra(:) ! at nit000 (Euler) -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
r8058 r8059 46 46 LOGICAL , PUBLIC :: ln_qsr_ice !: light penetration for ice-model LIM3 (clem) 47 47 INTEGER , PUBLIC :: nn_chldta !: use Chlorophyll data (=1) or not (=0) 48 INTEGER , PUBLIC :: nn_kd490dta !: use kd490dta data (=1) or not (=0) 48 49 REAL(wp), PUBLIC :: rn_abs !: fraction absorbed in the very near surface (RGB & 2 bands) 49 50 REAL(wp), PUBLIC :: rn_si0 !: very near surface depth of extinction (RGB & 2 bands) … … 54 55 REAL(wp) :: xsi1r !: inverse of rn_si1 55 56 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_chl ! structure of input Chl (file informations, fields read) 57 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_kd490 ! structure of input kd490 (file informations, fields read) 56 58 INTEGER, PUBLIC :: nksr ! levels below which the light cannot penetrate ( depth larger than 391 m) 57 59 REAL(wp), DIMENSION(3,61) :: rkrgb !: tabulated attenuation coefficients for RGB absorption … … 306 308 ! 307 309 ENDIF 310 ! slwa 311 IF( nn_kd490dta == 1 ) THEN ! use KD490 data read in ! 312 ! ! ------------------------- ! 313 nksr = jpk - 1 314 ! 315 CALL fld_read( kt, 1, sf_kd490 ) ! Read kd490 data and provide it at the current time step 316 ! 317 zcoef = ( 1. - rn_abs ) 318 ze0(:,:,1) = rn_abs * qsr(:,:) 319 ze1(:,:,1) = zcoef * qsr(:,:) 320 zea(:,:,1) = qsr(:,:) 321 ! 322 DO jk = 2, nksr+1 323 !CDIR NOVERRCHK 324 DO jj = 1, jpj 325 !CDIR NOVERRCHK 326 DO ji = 1, jpi 327 zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r ) 328 zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * sf_kd490(1)%fnow(ji,jj,1) ) 329 ze0(ji,jj,jk) = zc0 330 ze1(ji,jj,jk) = zc1 331 zea(ji,jj,jk) = ( zc0 + zc1 ) * tmask(ji,jj,jk) 332 END DO 333 END DO 334 END DO 335 ! clem: store attenuation coefficient of the first ocean level 336 IF ( ln_qsr_ice ) THEN 337 DO jj = 1, jpj 338 DO ji = 1, jpi 339 zzc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r ) 340 zzc1 = zcoef * EXP( - fse3t(ji,jj,1) * sf_kd490(1)%fnow(ji,jj,1) ) 341 fraqsr_1lev(ji,jj) = 1.0 - ( zzc0 + zzc1 ) * tmask(ji,jj,2) 342 END DO 343 END DO 344 ENDIF 345 ! 346 DO jk = 1, nksr ! compute and add qsr trend to ta 347 qsr_hc(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) 348 END DO 349 zea(:,:,nksr+1:jpk) = 0.e0 ! 350 CALL iom_put( 'qsr3d', zea ) ! Shortwave Radiation 3D distribution 351 ! 352 ENDIF ! use KD490 data 353 !slwa 308 354 ! 309 355 ! Add to the general trend … … 374 420 CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files 375 421 TYPE(FLD_N) :: sn_chl ! informations about the chlorofyl field to be read 376 !! 377 NAMELIST/namtra_qsr/ sn_chl, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice, & 378 & nn_chldta, rn_abs, rn_si0, rn_si1 422 TYPE(FLD_N) :: sn_kd490 ! informations about the kd490 field to be read 423 !! 424 NAMELIST/namtra_qsr/ sn_chl, sn_kd490, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice, & 425 & nn_chldta, rn_abs, rn_si0, rn_si1, nn_kd490dta 379 426 !!---------------------------------------------------------------------- 380 427 … … 409 456 WRITE(numout,*) ' RGB & 2 bands: shortess depth of extinction rn_si0 = ', rn_si0 410 457 WRITE(numout,*) ' 2 bands: longest depth of extinction rn_si1 = ', rn_si1 458 WRITE(numout,*) ' read in KD490 data nn_kd490dta = ', nn_kd490dta 411 459 ENDIF 412 460 … … 422 470 IF( ln_qsr_2bd ) ioptio = ioptio + 1 423 471 IF( ln_qsr_bio ) ioptio = ioptio + 1 472 IF( nn_kd490dta == 1 ) ioptio = ioptio + 1 424 473 ! 425 474 IF( ioptio /= 1 ) & … … 431 480 IF( ln_qsr_2bd ) nqsr = 3 432 481 IF( ln_qsr_bio ) nqsr = 4 482 IF( nn_kd490dta == 1 ) nqsr = 5 433 483 ! 434 484 IF(lwp) THEN ! Print the choice … … 438 488 IF( nqsr == 3 ) WRITE(numout,*) ' 2 bands light penetration' 439 489 IF( nqsr == 4 ) WRITE(numout,*) ' bio-model light penetration' 490 IF( nqsr == 5 ) WRITE(numout,*) ' KD490 light penetration' 440 491 ENDIF 441 492 ! … … 447 498 xsi0r = 1.e0 / rn_si0 448 499 xsi1r = 1.e0 / rn_si1 500 IF( nn_kd490dta == 1 ) THEN !* KD490 data : set sf_kd490 structure 501 IF(lwp) WRITE(numout,*) 502 IF(lwp) WRITE(numout,*) ' KD490 read in a file' 503 ALLOCATE( sf_kd490(1), STAT=ierror ) 504 IF( ierror > 0 ) THEN 505 CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_kd490 structure' ) ; RETURN 506 ENDIF 507 ALLOCATE( sf_kd490(1)%fnow(jpi,jpj,1) ) 508 IF( sn_kd490%ln_tint )ALLOCATE( sf_kd490(1)%fdta(jpi,jpj,1,2) ) 509 ! ! fill sf_kd490 with sn_kd490 and control print 510 CALL fld_fill( sf_kd490, (/ sn_kd490 /), cn_dir, 'tra_qsr_init', & 511 & 'Solar penetration function of read KD490', 'namtra_qsr' ) 449 512 ! ! ---------------------------------- ! 450 IF( ln_qsr_rgb ) THEN ! Red-Green-Blue light penetration !513 ELSEIF( ln_qsr_rgb ) THEN ! Red-Green-Blue light penetration ! 451 514 ! ! ---------------------------------- ! 452 515 ! -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
r8058 r8059 25 25 USE trd_oce ! trends: ocean variables 26 26 USE trdtra ! trends manager: tracers 27 USE tradwl ! solar radiation penetration (downwell method) 27 28 ! 28 29 USE in_out_manager ! I/O manager … … 138 139 139 140 !!gm IF( .NOT.ln_traqsr ) qsr(:,:) = 0.e0 ! no solar radiation penetration 140 IF( .NOT.ln_traqsr ) THEN ! no solar radiation penetration141 IF( .NOT.ln_traqsr .and. .NOT.ln_tradwl ) THEN ! no solar radiation penetration 141 142 qns(:,:) = qns(:,:) + qsr(:,:) ! total heat flux in qns 142 143 qsr(:,:) = 0.e0 ! qsr set to zero -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90
r8058 r8059 203 203 DO jj = 2, jpjm1 204 204 DO ji = fs_2, fs_jpim1 ! vector opt. 205 #if defined key_tracer_budget 206 ! ptrd(ji,jj,jk) = - ( pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik) ) * tmask(ji,jj,jk) 207 ptrd(ji,jj,jk) = - pf (ji,jj,jk) * tmask(ji,jj,jk) 208 #else 205 209 ptrd(ji,jj,jk) = - ( pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik) & 206 210 & - ( pun(ji,jj,jk) - pun(ji-ii,jj-ij,jk-ik) ) * ptn(ji,jj,jk) ) & 207 211 & / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) * tmask(ji,jj,jk) 212 #endif 208 213 END DO 209 214 END DO -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
r8058 r8059 18 18 USE phycst ! physical constants 19 19 USE iom ! I/O library 20 USE eosbn2 ! for zdf_mxl_zint 20 21 USE lib_mpp ! MPP library 21 22 USE wrk_nemo ! work arrays … … 27 28 28 29 PUBLIC zdf_mxl ! called by step.F90 29 PUBLIC zdf_mxl_alloc ! Used in zdf_tke_init30 30 31 31 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: nmln !: number of level in the mixed layer (used by TOP) … … 33 33 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlp !: mixed layer depth (rho=rho0+zdcrit) [m] 34 34 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlpt !: mixed layer depth at t-points [m] 35 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hmld_zint !: vertically-interpolated mixed layer depth [m] 36 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: htc_mld ! Heat content of hmld_zint 37 LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ll_found ! Is T_b to be found by interpolation ? 38 LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ll_belowml ! Flag points below mixed layer when ll_found=F 35 39 36 40 REAL(wp), PUBLIC :: rho_c = 0.01_wp !: density criterion for mixed layer depth 37 41 REAL(wp) :: avt_c = 5.e-4_wp ! Kz criterion for the turbocline depth 42 43 TYPE, PUBLIC :: MXL_ZINT !: Structure for MLD defs 44 INTEGER :: mld_type ! mixed layer type 45 REAL(wp) :: zref ! depth of initial T_ref 46 REAL(wp) :: dT_crit ! Critical temp diff 47 REAL(wp) :: iso_frac ! Fraction of rn_dT_crit used 48 END TYPE MXL_ZINT 38 49 39 50 !! * Substitutions … … 52 63 zdf_mxl_alloc = 0 ! set to zero if no array to be allocated 53 64 IF( .NOT. ALLOCATED( nmln ) ) THEN 54 ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc ) 65 ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj), & 66 & htc_mld(jpi,jpj), & 67 & ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc ) 55 68 ! 56 69 IF( lk_mpp ) CALL mpp_sum ( zdf_mxl_alloc ) … … 128 141 iikn = nmln(ji,jj) 129 142 imkt = mikt(ji,jj) 130 hmld (ji,jj) = ( fsdepw(ji,jj,iiki ) - fsdepw(ji,jj,imkt ) )* ssmask(ji,jj) ! Turbocline depth131 hmlp (ji,jj) = ( fsdepw(ji,jj,iikn ) - fsdepw(ji,jj, imkt) ) * ssmask(ji,jj) ! Mixed layer depth132 hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt ) )* ssmask(ji,jj) ! depth of the last T-point inside the mixed layer143 hmld (ji,jj) = ( fsdepw(ji,jj,iiki ) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj) ! Turbocline depth 144 hmlp (ji,jj) = ( fsdepw(ji,jj,iikn ) - fsdepw(ji,jj,MAX( imkt,nla10 ) ) ) * ssmask(ji,jj) ! Mixed layer depth 145 hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj) ! depth of the last T-point inside the mixed layer 133 146 END DO 134 147 END DO … … 138 151 ENDIF 139 152 153 ! Vertically-interpolated mixed-layer depth diagnostic 154 CALL zdf_mxl_zint( kt ) 155 140 156 IF(ln_ctl) CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 ) 141 157 ! … … 145 161 ! 146 162 END SUBROUTINE zdf_mxl 163 164 SUBROUTINE zdf_mxl_zint_mld( sf ) 165 !!---------------------------------------------------------------------------------- 166 !! *** ROUTINE zdf_mxl_zint_mld *** 167 ! 168 ! Calculate vertically-interpolated mixed layer depth diagnostic. 169 ! 170 ! This routine can calculate the mixed layer depth diagnostic suggested by 171 ! Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate 172 ! vertically-interpolated mixed-layer depth diagnostics with other parameter 173 ! settings set in the namzdf_mldzint namelist. 174 ! 175 ! If mld_type=1 the mixed layer depth is calculated as the depth at which the 176 ! density has increased by an amount equivalent to a temperature difference of 177 ! 0.8C at the surface. 178 ! 179 ! For other values of mld_type the mixed layer is calculated as the depth at 180 ! which the temperature differs by 0.8C from the surface temperature. 181 ! 182 ! David Acreman, Daley Calvert 183 ! 184 !!----------------------------------------------------------------------------------- 185 186 TYPE(MXL_ZINT), INTENT(in) :: sf 187 188 ! Diagnostic criteria 189 INTEGER :: nn_mld_type ! mixed layer type 190 REAL(wp) :: rn_zref ! depth of initial T_ref 191 REAL(wp) :: rn_dT_crit ! Critical temp diff 192 REAL(wp) :: rn_iso_frac ! Fraction of rn_dT_crit used 193 194 ! Local variables 195 REAL(wp), PARAMETER :: zepsilon = 1.e-30 ! local small value 196 INTEGER, POINTER, DIMENSION(:,:) :: ikmt ! number of active tracer levels 197 INTEGER, POINTER, DIMENSION(:,:) :: ik_ref ! index of reference level 198 INTEGER, POINTER, DIMENSION(:,:) :: ik_iso ! index of last uniform temp level 199 REAL, POINTER, DIMENSION(:,:,:) :: zT ! Temperature or density 200 REAL, POINTER, DIMENSION(:,:) :: ppzdep ! depth for use in calculating d(rho) 201 REAL, POINTER, DIMENSION(:,:) :: zT_ref ! reference temperature 202 REAL :: zT_b ! base temperature 203 REAL, POINTER, DIMENSION(:,:,:) :: zdTdz ! gradient of zT 204 REAL, POINTER, DIMENSION(:,:,:) :: zmoddT ! Absolute temperature difference 205 REAL :: zdz ! depth difference 206 REAL :: zdT ! temperature difference 207 REAL, POINTER, DIMENSION(:,:) :: zdelta_T ! difference critereon 208 REAL, POINTER, DIMENSION(:,:) :: zRHO1, zRHO2 ! Densities 209 INTEGER :: ji, jj, jk ! loop counter 210 211 !!------------------------------------------------------------------------------------- 212 ! 213 CALL wrk_alloc( jpi, jpj, ikmt, ik_ref, ik_iso) 214 CALL wrk_alloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 ) 215 CALL wrk_alloc( jpi, jpj, jpk, zT, zdTdz, zmoddT ) 216 217 ! Unpack structure 218 nn_mld_type = sf%mld_type 219 rn_zref = sf%zref 220 rn_dT_crit = sf%dT_crit 221 rn_iso_frac = sf%iso_frac 222 223 ! Set the mixed layer depth criterion at each grid point 224 IF( nn_mld_type == 0 ) THEN 225 zdelta_T(:,:) = rn_dT_crit 226 zT(:,:,:) = rhop(:,:,:) 227 ELSE IF( nn_mld_type == 1 ) THEN 228 ppzdep(:,:)=0.0 229 call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) ) 230 ! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST 231 ! [assumes number of tracers less than number of vertical levels] 232 zT(:,:,1:jpts)=tsn(:,:,1,1:jpts) 233 zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit 234 CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) ) 235 zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0 236 ! RHO from eos (2d version) doesn't calculate north or east halo: 237 CALL lbc_lnk( zdelta_T, 'T', 1. ) 238 zT(:,:,:) = rhop(:,:,:) 239 ELSE 240 zdelta_T(:,:) = rn_dT_crit 241 zT(:,:,:) = tsn(:,:,:,jp_tem) 242 END IF 243 244 ! Calculate the gradient of zT and absolute difference for use later 245 DO jk = 1 ,jpk-2 246 zdTdz(:,:,jk) = ( zT(:,:,jk+1) - zT(:,:,jk) ) / fse3w(:,:,jk+1) 247 zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) ) 248 END DO 249 250 ! Find density/temperature at the reference level (Kara et al use 10m). 251 ! ik_ref is the index of the box centre immediately above or at the reference level 252 ! Find rn_zref in the array of model level depths and find the ref 253 ! density/temperature by linear interpolation. 254 DO jk = jpkm1, 2, -1 255 WHERE ( fsdept(:,:,jk) > rn_zref ) 256 ik_ref(:,:) = jk - 1 257 zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - fsdept(:,:,jk-1) ) 258 END WHERE 259 END DO 260 261 ! If the first grid box centre is below the reference level then use the 262 ! top model level to get zT_ref 263 WHERE ( fsdept(:,:,1) > rn_zref ) 264 zT_ref = zT(:,:,1) 265 ik_ref = 1 266 END WHERE 267 268 ! The number of active tracer levels is 1 less than the number of active w levels 269 ikmt(:,:) = mbathy(:,:) - 1 270 271 ! Initialize / reset 272 ll_found(:,:) = .false. 273 274 IF ( rn_iso_frac - zepsilon > 0. ) THEN 275 ! Search for a uniform density/temperature region where adjacent levels 276 ! differ by less than rn_iso_frac * deltaT. 277 ! ik_iso is the index of the last level in the uniform layer 278 ! ll_found indicates whether the mixed layer depth can be found by interpolation 279 ik_iso(:,:) = ik_ref(:,:) 280 DO jj = 1, nlcj 281 DO ji = 1, nlci 282 !CDIR NOVECTOR 283 DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1 284 IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN 285 ik_iso(ji,jj) = jk 286 ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 287 EXIT 288 END IF 289 END DO 290 END DO 291 END DO 292 293 ! Use linear interpolation to find depth of mixed layer base where possible 294 hmld_zint(:,:) = rn_zref 295 DO jj = 1, jpj 296 DO ji = 1, jpi 297 IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN 298 zdz = abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 299 hmld_zint(ji,jj) = fsdept(ji,jj,ik_iso(ji,jj)) + zdz 300 END IF 301 END DO 302 END DO 303 END IF 304 305 ! If ll_found = .false. then calculate MLD using difference of zdelta_T 306 ! from the reference density/temperature 307 308 ! Prevent this section from working on land points 309 WHERE ( tmask(:,:,1) /= 1.0 ) 310 ll_found = .true. 311 END WHERE 312 313 DO jk=1, jpk 314 ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:) 315 END DO 316 317 ! Set default value where interpolation cannot be used (ll_found=false) 318 DO jj = 1, jpj 319 DO ji = 1, jpi 320 IF ( .not. ll_found(ji,jj) ) hmld_zint(ji,jj) = fsdept(ji,jj,ikmt(ji,jj)) 321 END DO 322 END DO 323 324 DO jj = 1, jpj 325 DO ji = 1, jpi 326 !CDIR NOVECTOR 327 DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj) 328 IF ( ll_found(ji,jj) ) EXIT 329 IF ( ll_belowml(ji,jj,jk) ) THEN 330 zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) ) 331 zdT = zT_b - zT(ji,jj,jk-1) 332 zdz = zdT / zdTdz(ji,jj,jk-1) 333 hmld_zint(ji,jj) = fsdept(ji,jj,jk-1) + zdz 334 EXIT 335 END IF 336 END DO 337 END DO 338 END DO 339 340 hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1) 341 ! 342 CALL wrk_dealloc( jpi, jpj, ikmt, ik_ref, ik_iso) 343 CALL wrk_dealloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 ) 344 CALL wrk_dealloc( jpi,jpj, jpk, zT, zdTdz, zmoddT ) 345 ! 346 END SUBROUTINE zdf_mxl_zint_mld 347 348 SUBROUTINE zdf_mxl_zint_htc( kt ) 349 !!---------------------------------------------------------------------- 350 !! *** ROUTINE zdf_mxl_zint_htc *** 351 !! 352 !! ** Purpose : 353 !! 354 !! ** Method : 355 !!---------------------------------------------------------------------- 356 357 INTEGER, INTENT(in) :: kt ! ocean time-step index 358 359 INTEGER :: ji, jj, jk 360 INTEGER :: ikmax 361 REAL(wp) :: zc, zcoef 362 ! 363 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ilevel 364 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zthick_0, zthick 365 366 !!---------------------------------------------------------------------- 367 368 IF( .NOT. ALLOCATED(ilevel) ) THEN 369 ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), & 370 & zthick(jpi,jpj), STAT=ji ) 371 IF( lk_mpp ) CALL mpp_sum(ji) 372 IF( ji /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' ) 373 ENDIF 374 375 ! Find last whole model T level above the MLD 376 ilevel(:,:) = 0 377 zthick_0(:,:) = 0._wp 378 379 DO jk = 1, jpkm1 380 DO jj = 1, jpj 381 DO ji = 1, jpi 382 zthick_0(ji,jj) = zthick_0(ji,jj) + fse3t(ji,jj,jk) 383 IF( zthick_0(ji,jj) < hmld_zint(ji,jj) ) ilevel(ji,jj) = jk 384 END DO 385 END DO 386 WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2) 387 WRITE(numout,*) 'fsdepw(jk+1 =',jk+1,') =',fsdepw(2,2,jk+1) 388 END DO 389 390 ! Surface boundary condition 391 IF( lk_vvl ) THEN ; zthick(:,:) = 0._wp ; htc_mld(:,:) = 0._wp 392 ELSE ; zthick(:,:) = sshn(:,:) ; htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1) 393 ENDIF 394 395 ! Deepest whole T level above the MLD 396 ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 ) 397 398 ! Integration down to last whole model T level 399 DO jk = 1, ikmax 400 DO jj = 1, jpj 401 DO ji = 1, jpi 402 zc = fse3t(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1 ) ) ! 0 below ilevel 403 zthick(ji,jj) = zthick(ji,jj) + zc 404 htc_mld(ji,jj) = htc_mld(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 405 END DO 406 END DO 407 END DO 408 409 ! Subsequent partial T level 410 zthick(:,:) = hmld_zint(:,:) - zthick(:,:) ! remaining thickness to reach MLD 411 412 DO jj = 1, jpj 413 DO ji = 1, jpi 414 htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem) & 415 & * MIN( fse3t(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1) 416 END DO 417 END DO 418 419 WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2) 420 421 ! Convert to heat content 422 zcoef = rau0 * rcp 423 htc_mld(:,:) = zcoef * htc_mld(:,:) 424 425 END SUBROUTINE zdf_mxl_zint_htc 426 427 SUBROUTINE zdf_mxl_zint( kt ) 428 !!---------------------------------------------------------------------- 429 !! *** ROUTINE zdf_mxl_zint *** 430 !! 431 !! ** Purpose : 432 !! 433 !! ** Method : 434 !!---------------------------------------------------------------------- 435 436 INTEGER, INTENT(in) :: kt ! ocean time-step index 437 438 INTEGER :: ios 439 INTEGER :: jn 440 441 INTEGER :: nn_mld_diag = 0 ! number of diagnostics 442 443 CHARACTER(len=1) :: cmld 444 445 TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 446 TYPE(MXL_ZINT), SAVE, DIMENSION(5) :: mld_diags 447 448 NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 449 450 !!---------------------------------------------------------------------- 451 452 IF( kt == nit000 ) THEN 453 REWIND( numnam_ref ) ! Namelist namzdf_mldzint in reference namelist 454 READ ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901) 455 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist', lwp ) 456 457 REWIND( numnam_cfg ) ! Namelist namzdf_mldzint in configuration namelist 458 READ ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 ) 459 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist', lwp ) 460 IF(lwm) WRITE ( numond, namzdf_mldzint ) 461 462 IF( nn_mld_diag > 5 ) CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' ) 463 464 mld_diags(1) = sn_mld1 465 mld_diags(2) = sn_mld2 466 mld_diags(3) = sn_mld3 467 mld_diags(4) = sn_mld4 468 mld_diags(5) = sn_mld5 469 470 IF( nn_mld_diag > 0 ) THEN 471 WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================' 472 WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)' 473 DO jn = 1, nn_mld_diag 474 WRITE(numout,*) 'MLD criterion',jn,':' 475 WRITE(numout,*) ' nn_mld_type =', mld_diags(jn)%mld_type 476 WRITE(numout,*) ' rn_zref =' , mld_diags(jn)%zref 477 WRITE(numout,*) ' rn_dT_crit =' , mld_diags(jn)%dT_crit 478 WRITE(numout,*) ' rn_iso_frac =', mld_diags(jn)%iso_frac 479 END DO 480 WRITE(numout,*) '====================================================================' 481 ENDIF 482 ENDIF 483 484 IF( nn_mld_diag > 0 ) THEN 485 DO jn = 1, nn_mld_diag 486 WRITE(cmld,'(I1)') jn 487 IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN 488 CALL zdf_mxl_zint_mld( mld_diags(jn) ) 489 490 IF( iom_use( "mldzint_"//cmld ) ) THEN 491 CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) ) 492 ENDIF 493 494 IF( iom_use( "mldhtc_"//cmld ) ) THEN 495 CALL zdf_mxl_zint_htc( kt ) 496 CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:) ) 497 ENDIF 498 ENDIF 499 END DO 500 ENDIF 501 502 END SUBROUTINE zdf_mxl_zint 147 503 148 504 !!====================================================================== -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r8058 r8059 85 85 USE stopar 86 86 USE stopts 87 USE diatmb ! Top,middle,bottom output 88 USE dia25h ! 25h mean output 87 89 88 90 IMPLICIT NONE … … 475 477 IF( lk_asminc ) CALL asm_inc_init ! Initialize assimilation increments 476 478 IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler 479 CALL dia_tmb_init ! TMB outputs 480 CALL dia_25h_init ! 25h mean outputs 477 481 ! 478 482 END SUBROUTINE nemo_init … … 630 634 USE ldftra_oce, ONLY: ldftra_oce_alloc 631 635 USE trc_oce , ONLY: trc_oce_alloc 636 USE diainsitutem, ONLY: insitu_tem_alloc 632 637 #if defined key_diadct 633 638 USE diadct , ONLY: diadct_alloc … … 646 651 ierr = ierr + ldftra_oce_alloc() ! ocean lateral physics : tracers 647 652 ierr = ierr + zdf_oce_alloc () ! ocean vertical physics 653 ierr = ierr + insitu_tem_alloc() 648 654 ! 649 655 ierr = ierr + trc_oce_alloc () ! shared TRC / TRA arrays -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/step.F90
r8058 r8059 255 255 CALL tra_sbc ( kstp ) ! surface boundary condition 256 256 IF( ln_traqsr ) CALL tra_qsr ( kstp ) ! penetrative solar radiation qsr 257 IF( ln_tradwl ) CALL tra_dwl ( kstp ) ! Polcoms Style Short Wave Radiation 257 258 IF( ln_trabbc ) CALL tra_bbc ( kstp ) ! bottom heat flux 258 259 IF( lk_trabbl ) CALL tra_bbl ( kstp ) ! advective (and/or diffusive) bottom boundary layer scheme … … 337 338 IF( lk_vvl ) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors 338 339 ! 339 IF( lrst_oce ) CALL rst_write( kstp ) ! write output ocean restart file340 340 341 341 #if defined key_agrif … … 361 361 CALL dia_wri_state( 'output.abort', kstp ) 362 362 ENDIF 363 IF( ln_harm_ana_store ) CALL harm_ana( kstp ) ! Harmonic analysis of tides 363 364 IF( kstp == nit000 ) THEN 364 365 CALL iom_close( numror ) ! close input ocean restart file … … 366 367 IF( lwm.AND.numoni /= -1 ) CALL FLUSH ( numoni ) ! flush output namelist ice 367 368 ENDIF 369 IF( lrst_oce ) CALL rst_write( kstp ) ! write output ocean restart file 368 370 369 371 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
r8058 r8059 25 25 USE sbcrnf ! surface boundary condition: runoff variables 26 26 USE sbccpl ! surface boundary condition: coupled formulation (call send at end of step) 27 USE sbcflx ! surface boundary condition: Fluxes 27 28 USE sbc_oce ! surface boundary condition: ocean 28 29 USE sbctide ! Tide initialisation … … 30 31 31 32 USE traqsr ! solar radiation penetration (tra_qsr routine) 33 USE tradwl ! POLCOMS style solar radiation (tra_dwl routine) 32 34 USE trasbc ! surface boundary condition (tra_sbc routine) 33 35 USE trabbc ! bottom boundary condition (tra_bbc routine) … … 106 108 USE prtctl ! Print control (prt_ctl routine) 107 109 110 USE harmonic_analysis ! harmonic analysis of tides (harm_ana routine) 111 USE bdytides ! harmonic analysis of tides (harm_ana routine) 108 112 USE diaobs ! Observation operator 109 113 -
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90
r8058 r8059 27 27 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: etot3 !: light absortion coefficient 28 28 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: facvol !: volume for degraded regions 29 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: rlambda2 !: Lambda2 for downwell version of Short wave Radiation 30 REAL(wp), PUBLIC :: rlambda !: Lambda for downwell version of Short wave Radiation 29 31 30 32 #if defined key_top … … 78 80 !! *** trc_oce_alloc *** 79 81 !!---------------------------------------------------------------------- 80 INTEGER :: ierr( 2) ! Local variables82 INTEGER :: ierr(3) ! Local variables 81 83 !!---------------------------------------------------------------------- 82 84 ierr(:) = 0 83 85 ALLOCATE( etot3 (jpi,jpj,jpk), STAT=ierr(1) ) 84 86 IF( lk_degrad) ALLOCATE( facvol(jpi,jpj,jpk), STAT=ierr(2) ) 87 ALLOCATE( rlambda2(jpi,jpj), STAT=ierr(3) ) 85 88 trc_oce_alloc = MAXVAL( ierr ) 86 89 ! 87 90 IF( trc_oce_alloc /= 0 ) CALL ctl_warn('trc_oce_alloc: failed to allocate etot3 array') 91 IF( trc_oce_alloc /= 0 ) CALL ctl_warn('trc_oce_alloc: failed to allocate etot3, facvol or rlambda2 array') 88 92 END FUNCTION trc_oce_alloc 89 93
Note: See TracChangeset
for help on using the changeset viewer.
