Changeset 6012 for branches/2015/dev_MetOffice_merge_2015/NEMOGCM
- Timestamp:
- 2015-12-07T16:11:45+01:00 (9 years ago)
- Location:
- branches/2015/dev_MetOffice_merge_2015/NEMOGCM
- Files:
-
- 33 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/SHARED/namelist_ref
r6008 r6012 245 245 ln_dm2dc = .false. ! daily mean to diurnal cycle on short wave 246 246 ln_rnf = .true. ! runoffs (T => fill namsbc_rnf) 247 nn_isf = 0 ! ice shelf melting/freezing (/=0 => fill namsbc_isf) 248 ! 0 =no isf 1 = presence of ISF 249 ! 2 = bg03 parametrisation 3 = rnf file for isf 250 ! 4 = ISF fwf specified 251 ! option 1 and 4 need ln_isfcav = .true. (domzgr) 247 ln_isf = .false. ! ice shelf (T => fill namsbc_isf) 252 248 ln_ssr = .true. ! Sea Surface Restoring on T and/or S (T => fill namsbc_ssr) 253 249 nn_fwb = 2 ! FreshWater Budget: =0 unchecked … … 431 427 &namsbc_isf ! Top boundary layer (ISF) 432 428 !----------------------------------------------------------------------- 433 ! ! file name ! frequency (hours) ! variable ! time interpol. ! clim ! 'yearly'/ ! weights ! rotation ! 434 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! 429 ! ! file name ! frequency (hours) ! variable ! time interpol. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! 430 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename ! 435 431 ! nn_isf == 4 436 sn_qisf = 'rnfisf' , -12 ,'sohflisf', .false. , .true. , 'yearly' , '' , '' 437 sn_fwfisf = 'rnfisf' , -12 ,'sowflisf', .false. , .true. , 'yearly' , '' , '' 432 sn_fwfisf = 'rnfisf' , -12 ,'sowflisf', .false. , .true. , 'yearly' , '' , '' , '' 438 433 ! nn_isf == 3 439 sn_rnfisf = 'runoffs' , -12 ,'sofwfisf', .false. , .true. , 'yearly' , '' ,''434 sn_rnfisf = 'rnfisf' , -12 ,'sofwfisf', .false. , .true. , 'yearly' , '' , '' , '' 440 435 ! nn_isf == 2 and 3 441 sn_depmax_isf = 'r unoffs' , -12 ,'sozisfmax' , .false. , .true. , 'yearly' , '' ,''442 sn_depmin_isf = 'r unoffs' , -12 ,'sozisfmin' , .false. , .true. , 'yearly' , '' ,''436 sn_depmax_isf = 'rnfisf' , -12 ,'sozisfmax' , .false. , .true. , 'yearly' , '' , '' , '' 437 sn_depmin_isf = 'rnfisf' , -12 ,'sozisfmin' , .false. , .true. , 'yearly' , '' , '' , '' 443 438 ! nn_isf == 2 444 sn_Leff_isf = 'rnfisf' , 0 ,'Leff' , .false. , .true. , 'yearly' , '' ,''439 sn_Leff_isf = 'rnfisf' , -12 ,'Leff' , .false. , .true. , 'yearly' , '' , '' , '' 445 440 ! for all case 446 ln_divisf = .true. ! apply isf melting as a mass flux or in the salinity trend. (maybe I should remove this option as for runoff?) 441 nn_isf = 1 ! ice shelf melting/freezing 442 ! 1 = presence of ISF 2 = bg03 parametrisation 443 ! 3 = rnf file for isf 4 = ISF fwf specified 444 ! option 1 and 4 need ln_isfcav = .true. (domzgr) 447 445 ! only for nn_isf = 1 or 2 448 446 rn_gammat0 = 1.0e-4 ! gammat coefficient used in blk formula 449 447 rn_gammas0 = 1.0e-4 ! gammas coefficient used in blk formula 448 ! only for nn_isf = 1 or 4 449 rn_hisf_tbl = 30. ! thickness of the top boundary layer (Losh et al. 2008) 450 ! 0 => thickness of the tbl = thickness of the first wet cell 450 451 ! only for nn_isf = 1 451 nn_isfblk = 1 ! 1 ISOMIP ; 2 conservative (3 equation formulation, Jenkins et al. 1991 ??) 452 rn_hisf_tbl = 30. ! thickness of the top boundary layer (Losh et al. 2008) 453 ! 0 => thickness of the tbl = thickness of the first wet cell 454 ln_conserve = .true. ! conservative case (take into account meltwater advection) 452 nn_isfblk = 1 ! 1 ISOMIP like: 2 equations formulation (Hunter et al., 2006) 453 ! 2 ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2015) 455 454 nn_gammablk = 1 ! 0 = cst Gammat (= gammat/s) 456 455 ! 1 = velocity dependend Gamma (u* * gammat/s) (Jenkins et al. 2010) 457 ! if you want to keep the cd as in global config, adjust rn_gammat0 to compensate 458 ! 2 = velocity and stability dependent Gamma Holland et al. 1999 456 ! 2 = velocity and stability dependent Gamma (Holland et al. 1999) 459 457 / 460 458 !----------------------------------------------------------------------- -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/cfg.txt
r6006 r6012 7 7 AMM12 OPA_SRC 8 8 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 9 GYRE OPA_SRC10 9 ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC 11 10 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 12 11 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC 12 GYRE OPA_SRC -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/uspcfg.txt
r5333 r6012 1 1 ORCA1_CICE # ORCA2_LIM # OPA_SRC TOP_SRC # http://gws-access.ceda.ac.uk/public/nemo/uspconfigs/ORCA1_CICE/v3.6.0/ORCA1_CICE_ctl.txt 2 2 ISOMIP # GYRE # OPA_SRC # http://gws-access.ceda.ac.uk/public/nemo/uspconfigs/ISOMIP/v3.6.0/ISOMIP_ctl.txt 3 ISOMIP_r5151_UKMO_ISF # GYRE # OPA_SRC # http://gws-access.ceda.ac.uk/public/nemo/uspconfigs/ISOMIP_r5151_UKMO_ISF/v3.6.0/ISOMIP_ctl.txt -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90
r5930 r6012 489 489 CALL bn2 ( pts, rab_n, rn2 ) ! now Brunt-Vaisala 490 490 491 ! Partial steps: before Horizontal DErivative 492 IF( ln_zps .AND. .NOT. ln_isfcav) & 493 & CALL zps_hde ( kt, jpts, pts, gtsu, gtsv, & ! Partial steps: before horizontal gradient 494 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 495 IF( ln_zps .AND. ln_isfcav) & 496 & CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, & ! Partial steps for top cell (ISF) 497 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 498 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the first ocean level 491 ! Partial steps: before Horizontal DErivative 492 IF( ln_zps .AND. .NOT. ln_isfcav) & 493 & CALL zps_hde ( kt, jpts, pts, gtsu, gtsv, & ! Partial steps: before horizontal gradient 494 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 495 IF( ln_zps .AND. ln_isfcav) & 496 & CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 497 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level 499 498 500 499 rn2b(:,:,:) = rn2(:,:,:) ! need for zdfmxl -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r5997 r6012 636 636 !!gm 637 637 638 IF( ln_zps .AND. .NOT. lk_c1d ) THEN ! Partial steps: before horizontal gradient 639 IF(ln_isfcav) THEN ! ocean cavities: top and bottom cells (ISF) 640 CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi, & 641 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 642 & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) 643 ELSE ! no ocean cavities: bottom cells 644 CALL zps_hde ( kt, jpts, tsb, gtsu, gtsv, & ! 645 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 646 ENDIF 647 ENDIF 648 ! 638 IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav) & 639 & CALL zps_hde ( kt, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 640 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 641 IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) & 642 & CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 643 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the last ocean level 644 645 #if defined key_zdfkpp 646 CALL eos( tsn, rhd, fsdept_n(:,:,:) ) ! Compute rhd 647 !!gm fabien CALL eos( tsn, rhd ) ! Compute rhd 648 #endif 649 649 650 DEALLOCATE( t_bkginc ) 650 651 DEALLOCATE( s_bkginc ) -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90
r5930 r6012 93 93 ! ----------------------------------------------------------------------- 94 94 !!gm replace these lines : 95 z_cflxemp = SUM ( ( emp(:,:) -rnf(:,:)+fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rau095 z_cflxemp = SUM ( ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rau0 96 96 IF( lk_mpp ) CALL mpp_sum( z_cflxemp ) ! sum over the global domain 97 97 !!gm by : -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
r6006 r6012 194 194 DO ji = 1,jpi 195 195 ! Elevation 196 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)* tmask_i(ji,jj)196 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj) 197 197 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj) 198 198 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj) … … 358 358 X1=ana_amp(ji,jj,jh,1) 359 359 X2=-ana_amp(ji,jj,jh,2) 360 out_v(ji,jj, jh)=X1 * vmask_i(ji,jj)361 out_v(ji,jj,nb_ana+jh)=X2 * vmask_i(ji,jj)360 out_v(ji,jj, jh)=X1 * ssvmask(ji,jj) 361 out_v(ji,jj,nb_ana+jh)=X2 * ssvmask(ji,jj) 362 362 END DO 363 363 END DO -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
r6006 r6012 101 101 IF( ln_rnf_sal) z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * surf(:,:) ) 102 102 ! Add ice shelf heat & salt input 103 IF( nn_isf .GE. 1 ) THEN 104 z_frc_trd_t = z_frc_trd_t + glob_sum( risf_tsc(:,:,jp_tem) * surf(:,:) ) 105 z_frc_trd_s = z_frc_trd_s + glob_sum( risf_tsc(:,:,jp_sal) * surf(:,:) ) 106 ENDIF 107 103 IF( ln_isf ) z_frc_trd_t = z_frc_trd_t + glob_sum( risf_tsc(:,:,jp_tem) * surf(:,:) ) 108 104 ! Add penetrative solar radiation 109 105 IF( ln_traqsr ) z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr (:,:) * surf(:,:) ) -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r6006 r6012 256 256 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i !: interior domain T-point mask 257 257 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_h !: internal domain T-point mask (Figure 8.5 NEMO book) 258 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bmask 258 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bmask !: land/ocean mask of barotropic stream function 259 259 260 260 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfdep !: top first ocean level (ISF) … … 390 390 391 391 ALLOCATE( mbathy(jpi,jpj) , bathy (jpi,jpj) , & 392 & tmask_i(jpi,jpj) , tmask_h(jpi, jpj), &392 & tmask_i(jpi,jpj) , tmask_h(jpi,jpj) , & 393 393 & ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , & 394 394 & bmask(jpi,jpj) , & -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r5836 r6012 506 506 CALL iom_close( inum ) 507 507 508 !!gm THIS is TO BE REMOVED !!!!!!!509 510 ! need to be define for the extended grid south of -80S511 ! some point are undefined but you need to have e1 and e2 .NE. 0512 WHERE (e1t==0.0_wp)513 e1t=1.0e2514 END WHERE515 WHERE (e1v==0.0_wp)516 e1v=1.0e2517 END WHERE518 WHERE (e1u==0.0_wp)519 e1u=1.0e2520 END WHERE521 WHERE (e1f==0.0_wp)522 e1f=1.0e2523 END WHERE524 WHERE (e2t==0.0_wp)525 e2t=1.0e2526 END WHERE527 WHERE (e2v==0.0_wp)528 e2v=1.0e2529 END WHERE530 WHERE (e2u==0.0_wp)531 e2u=1.0e2532 END WHERE533 WHERE (e2f==0.0_wp)534 e2f=1.0e2535 END WHERE536 !!gm end537 538 508 END SUBROUTINE hgr_read 539 509 -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r6006 r6012 975 975 ! 976 976 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 977 IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )977 IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 978 978 ! 979 979 IF(lwp) THEN ! Print the choice -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r6006 r6012 1258 1258 END WHERE 1259 1259 1260 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1261 WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1) 1262 misfdep = 0; risfdep = 0.0_wp; 1263 mbathy = 0; bathy = 0.0_wp; 1264 END WHERE 1265 WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp) 1266 misfdep = 0; risfdep = 0.0_wp; 1267 mbathy = 0; bathy = 0.0_wp; 1268 END WHERE 1269 1260 1270 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation 1261 1271 icompt = 0 … … 1323 1333 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1324 1334 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1325 1326 1335 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1327 1336 IF (zbathydiff <= zrisfdepdiff) THEN … … 1369 1378 END DO 1370 1379 1371 ! point V mbathy(ji,jj) EQmisfdep(ji,jj+1)1380 ! point V mbathy(ji,jj) == misfdep(ji,jj+1) 1372 1381 DO jj = 1, jpjm1 1373 1382 DO ji = 1, jpim1 … … 1405 1414 mbathy(:,:) = INT( zbathy(:,:) ) 1406 1415 ENDIF 1407 ! point V misdep(ji,jj) EQmbathy(ji,jj+1)1416 ! point V misdep(ji,jj) == mbathy(ji,jj+1) 1408 1417 DO jj = 1, jpjm1 1409 1418 DO ji = 1, jpim1 … … 1443 1452 ENDIF 1444 1453 1445 ! point U mbathy(ji,jj) EQmisfdep(ji,jj+1)1454 ! point U mbathy(ji,jj) == misfdep(ji,jj+1) 1446 1455 DO jj = 1, jpjm1 1447 1456 DO ji = 1, jpim1 … … 1480 1489 ENDIF 1481 1490 1482 ! point U misfdep(ji,jj) EQbathy(ji,jj+1)1491 ! point U misfdep(ji,jj) == bathy(ji,jj+1) 1483 1492 DO jj = 1, jpjm1 1484 1493 DO ji = 1, jpim1 … … 1749 1758 ! end check compatibility ice shelf/bathy 1750 1759 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1760 WHERE (risfdep(:,:) <= 10._wp) 1761 misfdep = 1; risfdep = 0.0_wp; 1762 END WHERE 1763 1751 1764 IF( icompt == 0 ) THEN 1752 1765 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry' … … 1832 1845 ENDIF 1833 1846 ! ... on ik / ik-1 1834 e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik) 1847 e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1835 1848 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1836 1849 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90
r6006 r6012 20 20 USE oce ! ocean dynamics and tracers 21 21 USE dom_oce ! ocean space and time domain 22 USE sbc_oce, ONLY : ln_rnf, nn_isf ! surface boundary condition: ocean22 USE sbc_oce, ONLY : ln_rnf, ln_isf ! surface boundary condition: ocean 23 23 USE sbcrnf ! river runoff 24 24 USE sbcisf ! ice shelf … … 91 91 END DO 92 92 ! 93 IF( ln_rnf 93 IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) !== runoffs ==! (update hdivn field) 94 94 ! 95 IF( ln_ divisf .AND. nn_isf > 0) CALL sbc_isf_div( hdivn ) !== ice shelf ==! (update hdivn field)95 IF( ln_isf ) CALL sbc_isf_div( hdivn ) !== ice shelf ==! (update hdivn field) 96 96 ! 97 IF( ln_iscpl .AND. ln_hsb ) CALL iscpl_div( hdivn )!== ice sheet ==! (update hdivn field)97 IF( ln_iscpl .AND. ln_hsb ) CALL iscpl_div( hdivn ) !== ice sheet ==! (update hdivn field) 98 98 ! 99 CALL lbc_lnk( hdivn, 'T', 1. ) 99 CALL lbc_lnk( hdivn, 'T', 1. ) !== lateral boundary cond. ==! (no sign change) 100 100 ! 101 101 IF( nn_timing == 1 ) CALL timing_stop('div_hor') -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5930 r6012 45 45 USE wrk_nemo ! Memory Allocation 46 46 USE timing ! Timing 47 USE iom 47 48 48 49 IMPLICIT NONE … … 130 131 INTEGER :: ioptio = 0 ! temporary integer 131 132 INTEGER :: ios ! Local integer output status for namelist read 133 !! 134 INTEGER :: ji, jj, jk, ikt ! dummy loop indices ISF 135 REAL(wp) :: znad 136 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop, zrhd ! hypothesys on isf density 137 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_isf ! density at bottom of ISF 138 REAL(wp), POINTER, DIMENSION(:,:) :: ziceload ! density at bottom of ISF 132 139 !! 133 140 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & … … 190 197 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 191 198 ! 192 ! initialisation of ice load 193 riceload(:,:)=0.0 199 ! initialisation of ice shelf load 200 IF ( .NOT. ln_isfcav ) riceload(:,:)=0.0 201 IF ( ln_isfcav ) THEN 202 CALL wrk_alloc( jpi,jpj, 2, ztstop) 203 CALL wrk_alloc( jpi,jpj,jpk, zrhd ) 204 CALL wrk_alloc( jpi,jpj, zrhdtop_isf, ziceload) 205 ! 206 IF(lwp) WRITE(numout,*) 207 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 208 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 209 210 ! To use density and not density anomaly 211 znad=1._wp 212 213 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 214 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 215 216 ! compute density of the water displaced by the ice shelf 217 DO jk = 1, jpk 218 CALL eos(ztstop(:,:,:),fsdept_n(:,:,jk),zrhd(:,:,jk)) 219 END DO 220 221 ! compute rhd at the ice/oce interface (ice shelf side) 222 CALL eos(ztstop,risfdep,zrhdtop_isf) 223 224 ! Surface value + ice shelf gradient 225 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 226 ! divided by 2 later 227 ziceload = 0._wp 228 DO jj = 1, jpj 229 DO ji = 1, jpi 230 ikt=mikt(ji,jj) 231 ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 232 DO jk=2,ikt-1 233 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 234 & * (1._wp - tmask(ji,jj,jk)) 235 END DO 236 IF (ikt >= 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 237 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 238 END DO 239 END DO 240 riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5 241 242 CALL wrk_dealloc( jpi,jpj, 2, ztstop) 243 CALL wrk_dealloc( jpi,jpj,jpk, zrhd ) 244 CALL wrk_dealloc( jpi,jpj, zrhdtop_isf, ziceload) 245 END IF 194 246 ! 195 247 END SUBROUTINE dyn_hpg_init … … 447 499 SUBROUTINE hpg_isf( kt ) 448 500 !!--------------------------------------------------------------------- 449 !! *** ROUTINE hpg_ sco***501 !! *** ROUTINE hpg_isf *** 450 502 !! 451 503 !! ** Method : s-coordinate case. Jacobian scheme. … … 466 518 INTEGER, INTENT(in) :: kt ! ocean time-step index 467 519 !! 468 INTEGER :: ji, jj, jk, ik u, ikv, ikt, iktp1i, iktp1j! dummy loop indices469 REAL(wp) :: zcoef0, zuap, zvap, znad , ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1! temporary scalars470 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj , zrhd520 INTEGER :: ji, jj, jk, ikt, iktp1i, iktp1j ! dummy loop indices 521 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 522 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 471 523 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop 472 REAL(wp), POINTER, DIMENSION(:,:) :: ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj 473 !!---------------------------------------------------------------------- 474 ! 475 CALL wrk_alloc( jpi,jpj, 2, ztstop) 476 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 477 CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 478 ! 479 IF( kt == nit000 ) THEN 480 IF(lwp) WRITE(numout,*) 481 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 482 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 483 ENDIF 484 524 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_oce 525 !!---------------------------------------------------------------------- 526 ! 527 CALL wrk_alloc( jpi,jpj, 2, ztstop) 528 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj) 529 CALL wrk_alloc( jpi,jpj, zrhdtop_oce ) 530 ! 485 531 ! Local constant initialization 486 532 zcoef0 = - grav * 0.5_wp 533 487 534 ! To use density and not density anomaly 488 ! IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume489 ! ELSE ; znad = 0._wp ! Fixed volume490 ! ENDIF491 535 znad=1._wp 536 492 537 ! iniitialised to 0. zhpi zhpi 493 538 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 494 539 495 ! Partial steps: top & bottom before horizontal gradient496 !jc CALL zps_hde_isf( kt, jpts, tsn, gtsu, gtsv, gtui, gtvi, &497 !jc & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , &498 !jc & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi )499 500 !==================================================================================501 !=====Compute iceload and contribution of the half first wet layer =================502 !===================================================================================503 504 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)505 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp506 507 ! compute density of the water displaced by the ice shelf508 zrhd = rhd ! save rhd509 DO jk = 1, jpk510 zdept(:,:)=gdept_1d(jk)511 CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk))512 END DO513 WHERE ( tmask(:,:,:) == 1._wp)514 rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd515 END WHERE516 517 ! compute rhd at the ice/oce interface (ice shelf side)518 CALL eos(ztstop,risfdep,zrhdtop_isf)519 520 540 ! compute rhd at the ice/oce interface (ocean side) 541 ! usefull to reduce residual current in the test case ISOMIP with no melting 521 542 DO ji=1,jpi 522 543 DO jj=1,jpj … … 527 548 END DO 528 549 CALL eos(ztstop,risfdep,zrhdtop_oce) 529 ! 530 ! Surface value + ice shelf gradient 531 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 532 ziceload = 0._wp 533 DO jj = 1, jpj 534 DO ji = 1, jpi ! vector opt. 535 ikt=mikt(ji,jj) 536 ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 537 DO jk=2,ikt-1 538 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 539 & * (1._wp - tmask(ji,jj,jk)) 540 END DO 541 IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) & 542 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 543 END DO 544 END DO 545 riceload(:,:) = 0.0_wp ; riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5 546 ! compute zp from z=0 to first T wet point (correction due to zps not yet applied) 550 551 !================================================================================== 552 !===== Compute surface value ===================================================== 553 !================================================================================== 547 554 DO jj = 2, jpjm1 548 555 DO ji = fs_2, fs_jpim1 ! vector opt. … … 550 557 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 551 558 ! we assume ISF is in isostatic equilibrium 552 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj 553 & * ( 2._wp * znad + rhd(ji+1,jj ,iktp1i) + zrhdtop_oce(ji+1,jj) ) &554 & - 0.5_wp * fse3w(ji ,jj ,ikt )&555 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) )&556 & + ( ziceload(ji+1,jj) - ziceload(ji,jj)))557 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji 558 & * ( 2._wp * znad + rhd(ji ,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) &559 & - 0.5_wp * fse3w(ji ,jj ,ikt )&560 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) )&561 & + ( ziceload(ji,jj+1) - ziceload(ji,jj) ))559 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj,iktp1i) & 560 & * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) ) & 561 & - 0.5_wp * fse3w(ji,jj,ikt) & 562 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 563 & + ( riceload(ji+1,jj) - riceload(ji,jj)) ) 564 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji,jj+1,iktp1j) & 565 & * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 566 & - 0.5_wp * fse3w(ji,jj,ikt) & 567 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 568 & + ( riceload(ji,jj+1) - riceload(ji,jj)) ) 562 569 ! s-coordinate pressure gradient correction (=0 if z coordinate) 563 570 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & … … 570 577 END DO 571 578 END DO 572 !==================================================================================573 !===== Compute partial cell contribution for the top cell =========================574 !==================================================================================575 DO jj = 2, jpjm1576 DO ji = fs_2, fs_jpim1 ! vector opt.577 iku = miku(ji,jj) ;578 zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp579 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))580 ! u direction581 IF ( iku .GT. 1 ) THEN582 ! case iku583 zhpi(ji,jj,iku) = zcoef0 / e1u(ji,jj) * ze3wu &584 & * ( rhd (ji+1,jj,iku) + rhd (ji,jj,iku) &585 & + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad )586 ! corrective term ( = 0 if z coordinate )587 zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj)588 ! zhpi will be added in interior loop589 ua(ji,jj,iku) = ua(ji,jj,iku) + zuap590 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure591 IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj) = zhpi(ji,jj,iku)592 593 ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps)594 zhpiint = zcoef0 / e1u(ji,jj) &595 & * ( fse3w(ji+1,jj ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad) &596 & + (rhd(ji+1,jj,iku ) + znad) ) * tmask(ji+1,jj,iku) &597 & - fse3w(ji ,jj ,iku+1) * ( (rhd(ji ,jj,iku+1) + znad) &598 & + (rhd(ji ,jj,iku ) + znad) ) * tmask(ji ,jj,iku) )599 zhpi(ji,jj,iku+1) = zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint600 END IF601 602 ! v direction603 ikv = mikv(ji,jj)604 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))605 IF ( ikv .GT. 1 ) THEN606 ! case ikv607 zhpj(ji,jj,ikv) = zcoef0 / e2v(ji,jj) * ze3wv &608 & * ( rhd(ji,jj+1,ikv) + rhd (ji,jj,ikv) &609 & + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad )610 ! corrective term ( = 0 if z coordinate )611 zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj)612 ! zhpi will be added in interior loop613 va(ji,jj,ikv) = va(ji,jj,ikv) + zvap614 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure615 IF (mbkv(ji,jj) == ikv + 1) zpshpj(ji,jj) = zhpj(ji,jj,ikv)616 617 ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps)618 zhpjint = zcoef0 / e2v(ji,jj) &619 & * ( fse3w(ji ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad) &620 & + (rhd(ji,jj+1,ikv ) + znad) ) * tmask(ji,jj+1,ikv) &621 & - fse3w(ji ,jj ,ikv+1) * ( (rhd(ji,jj ,ikv+1) + znad) &622 & + (rhd(ji,jj ,ikv ) + znad) ) * tmask(ji,jj ,ikv) )623 zhpj(ji,jj,ikv+1) = zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint624 END IF625 END DO626 END DO627 579 628 580 !================================================================================== 629 581 !===== Compute interior value ===================================================== 630 582 !================================================================================== 631 632 DO jj = 2, jpjm1 633 DO ji = fs_2, fs_jpim1 ! vector opt. 634 iku=miku(ji,jj); ikv=mikv(ji,jj) 635 DO jk = 2, jpkm1 583 ! interior value (2=<jk=<jpkm1) 584 DO jk = 2, jpkm1 585 DO jj = 2, jpjm1 586 DO ji = fs_2, fs_jpim1 ! vector opt. 636 587 ! hydrostatic pressure gradient along s-surfaces 637 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 638 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1) & 639 & + zcoef0 / e1u(ji,jj) & 640 & * ( fse3w(ji+1,jj ,jk) * ( (rhd(ji+1,jj,jk ) + znad) & 641 & + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1) & 642 & - fse3w(ji ,jj ,jk) * ( (rhd(ji ,jj,jk ) + znad) & 643 & + (rhd(ji ,jj,jk-1) + znad) ) * tmask(ji ,jj,jk-1) ) 588 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 589 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 590 & - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 591 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 592 & * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 593 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 644 594 ! s-coordinate pressure gradient correction 645 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 646 zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 647 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1) 648 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 649 650 ! hydrostatic pressure gradient along s-surfaces 651 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 652 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1) & 653 & + zcoef0 / e2v(ji,jj) & 654 & * ( fse3w(ji ,jj+1,jk) * ( (rhd(ji,jj+1,jk ) + znad) & 655 & + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1) & 656 & - fse3w(ji ,jj ,jk) * ( (rhd(ji,jj ,jk ) + znad) & 657 & + (rhd(ji,jj ,jk-1) + znad) ) * tmask(ji,jj ,jk-1) ) 658 ! s-coordinate pressure gradient correction 659 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 660 zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 661 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1) 595 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 596 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 597 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 598 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 662 599 ! add to the general momentum trend 663 va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 600 ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 601 va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 664 602 END DO 665 603 END DO 666 604 END DO 667 668 !================================================================================== 669 !===== Compute bottom cell contribution (partial cell) ============================ 670 !================================================================================== 671 672 DO jj = 2, jpjm1 673 DO ji = 2, jpim1 674 iku = mbku(ji,jj) 675 ikv = mbkv(ji,jj) 676 677 IF (iku .GT. 1) THEN 678 ! remove old value (interior case) 679 zuap = -zcoef0 * ( rhd (ji+1,jj ,iku) + rhd (ji,jj,iku) + 2._wp * znad ) & 680 & * ( fsde3w(ji+1,jj ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj) 681 ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 682 ! put new value 683 ! -zpshpi to avoid double contribution of the partial step in the top layer 684 zuap = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj) / e1u(ji,jj) 685 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 686 ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 687 END IF 688 ! v direction 689 IF (ikv .GT. 1) THEN 690 ! remove old value (interior case) 691 zvap = -zcoef0 * ( rhd (ji ,jj+1,ikv) + rhd (ji,jj,ikv) + 2._wp * znad ) & 692 & * ( fsde3w(ji ,jj+1,ikv) - fsde3w(ji,jj,ikv) ) / e2v(ji,jj) 693 va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 694 ! put new value 695 ! -zpshpj to avoid double contribution of the partial step in the top layer 696 zvap = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) / e2v(ji,jj) 697 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj) 698 va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 699 END IF 700 END DO 701 END DO 702 703 ! set back to original density value into the ice shelf cell (maybe useless because it is masked) 704 rhd = zrhd 705 ! 706 CALL wrk_dealloc( jpi,jpj,2, ztstop) 707 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 708 CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 605 ! 606 CALL wrk_dealloc( jpi,jpj,2 , ztstop) 607 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj) 608 CALL wrk_dealloc( jpi,jpj , zrhdtop_oce ) 709 609 ! 710 610 END SUBROUTINE hpg_isf -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r6006 r6012 95 95 ! 96 96 INTEGER :: ji, jj, jk ! dummy loop indices 97 INTEGER :: ikt ! local integers 97 98 REAL(wp) :: zue3a, zue3n, zue3b, zuf ! local scalars 98 99 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - … … 220 221 ! Add volume filter correction: compatibility with tracer advection scheme 221 222 ! => time filter + conservation correction (only at the first level) 222 IF ( nn_isf == 0) THEN ! if no ice shelf melting223 IF ( .NOT. ln_isf ) THEN ! if no ice shelf melting 223 224 fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) & 224 & -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)225 & - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 225 226 ELSE ! if ice shelf melting 226 227 DO jj = 1,jpj 227 228 DO ji = 1,jpi 228 jk= mikt(ji,jj)229 fse3t_b(ji,jj, jk) = fse3t_b(ji,jj,jk) - atfp * rdt * r1_rau0&229 ikt = mikt(ji,jj) 230 fse3t_b(ji,jj,ikt) = fse3t_b(ji,jj,ikt) - atfp * rdt * r1_rau0 & 230 231 & * ( (emp_b(ji,jj) - emp(ji,jj) ) & 231 232 & - (rnf_b(ji,jj) - rnf(ji,jj) ) & 232 & + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj, jk)233 & + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,ikt) 233 234 END DO 234 235 END DO -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r5930 r6012 231 231 IF( ioptio > 1 .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) ) & 232 232 & CALL ctl_stop( ' Choose only one surface pressure gradient scheme' ) 233 IF( ln_dynspg_ ts.AND. ln_isfcav ) &234 & CALL ctl_stop( ' dynspg_ tsnot tested with ice shelf cavity ' )233 IF( ln_dynspg_exp .AND. ln_isfcav ) & 234 & CALL ctl_stop( ' dynspg_exp not tested with ice shelf cavity ' ) 235 235 ! 236 236 IF( ln_dynspg_exp) nspg = 0 -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5930 r6012 145 145 INTEGER :: ji, jj, jk, jn ! dummy loop indices 146 146 INTEGER :: ikbu, ikbv, noffset ! local integers 147 INTEGER :: iktu, iktv ! local integers 147 148 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 148 149 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - … … 384 385 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 385 386 DO ji = fs_2, fs_jpim1 386 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)387 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)387 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 388 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 388 389 END DO 389 390 END DO … … 413 414 zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 414 415 zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 416 ! 417 ! ! Add top stress contribution from baroclinic velocities: 418 IF (ln_bt_fw) THEN 419 DO jj = 2, jpjm1 420 DO ji = fs_2, fs_jpim1 ! vector opt. 421 iktu = miku(ji,jj) 422 iktv = mikv(ji,jj) 423 zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities 424 zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 425 END DO 426 END DO 427 ELSE 428 DO jj = 2, jpjm1 429 DO ji = fs_2, fs_jpim1 ! vector opt. 430 iktu = miku(ji,jj) 431 iktv = mikv(ji,jj) 432 zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities 433 zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 434 END DO 435 END DO 436 ENDIF 437 ! 438 ! Note that the "unclipped" top friction parameter is used even with explicit drag 439 zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * tfrua(:,:) * zwx(:,:) 440 zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * tfrva(:,:) * zwy(:,:) 415 441 ! 416 442 IF (ln_bt_fw) THEN ! Add wind forcing … … 544 570 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 545 571 DO ji = 2, fs_jpim1 ! Vector opt. 546 zwx(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) &572 zwx(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 547 573 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 548 574 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 549 zwy(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) &575 zwy(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 550 576 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 551 577 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) … … 607 633 END DO 608 634 END DO 609 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1)635 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 610 636 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 611 637 … … 622 648 DO jj = 2, jpjm1 623 649 DO ji = 2, jpim1 ! NO Vector Opt. 624 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj)&625 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &626 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )627 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj)&628 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &629 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )650 zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 651 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 652 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 653 zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 654 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 655 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 630 656 END DO 631 657 END DO … … 661 687 DO jj = 2, jpjm1 662 688 DO ji = 2, jpim1 663 zx1 = z1_2 * umask(ji ,jj,1) * r1_e1e2u(ji ,jj) &689 zx1 = z1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) & 664 690 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) & 665 691 & + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) ) 666 zy1 = z1_2 * vmask(ji ,jj,1) * r1_e1e2v(ji ,jj ) &692 zy1 = z1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) & 667 693 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) & 668 694 & + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) … … 736 762 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 737 763 ! 764 ! Add top stresses: 765 zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:) 766 zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 767 ! 738 768 ! Surface pressure trend: 739 769 DO jj = 2, jpjm1 … … 755 785 & + zu_trd(ji,jj) & 756 786 & + zu_frc(ji,jj) ) & 757 & ) * umask(ji,jj,1)787 & ) * ssumask(ji,jj) 758 788 759 789 va_e(ji,jj) = ( vn_e(ji,jj) & … … 761 791 & + zv_trd(ji,jj) & 762 792 & + zv_frc(ji,jj) ) & 763 & ) * vmask(ji,jj,1)764 END DO 765 END DO 766 767 ELSE 793 & ) * ssvmask(ji,jj) 794 END DO 795 END DO 796 797 ELSE ! Flux form 768 798 DO jj = 2, jpjm1 769 799 DO ji = fs_2, fs_jpim1 ! vector opt. 770 800 771 zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1))772 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))801 zhura = ssumask(ji,jj)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj)) 802 zhvra = ssvmask(ji,jj)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj)) 773 803 774 804 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & … … 791 821 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 792 822 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 793 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )794 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )823 hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 824 hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) 795 825 ! 796 826 ENDIF … … 827 857 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 828 858 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 829 ELSE 859 ELSE ! Sum transports 830 860 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 831 861 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) … … 881 911 END DO 882 912 ! Save barotropic velocities not transport: 883 ua_b (:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )884 va_b (:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )913 ua_b (:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 914 va_b (:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 885 915 ENDIF 886 916 ! … … 897 927 ! 898 928 IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN 899 IF ( Agrif_NbStepint() .EQ.0 ) THEN929 IF ( Agrif_NbStepint() == 0 ) THEN 900 930 ub2_i_b(:,:) = 0.e0 901 931 vb2_i_b(:,:) = 0.e0 -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
r5836 r6012 90 90 DO jj = 2, jpjm1 ! vertical momentum advection at w-point 91 91 DO ji = fs_2, fs_jpim1 ! vector opt. 92 zwuw(ji,jj,jk) = ( zww(ji+1,jj ) + zww(ji,jj) ) * ( un(ji,jj,jk-1) -un(ji,jj,jk) )93 zwvw(ji,jj,jk) = ( zww(ji ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1) -vn(ji,jj,jk) )92 zwuw(ji,jj,jk) = ( zww(ji+1,jj ) + zww(ji,jj) ) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) 93 zwvw(ji,jj,jk) = ( zww(ji ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) 94 94 END DO 95 95 END DO -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r5836 r6012 112 112 !! 113 113 INTEGER :: ji , jj , jk ! dummy loop indices 114 INTEGER :: ii0, ii1 , iku! temporary integer115 INTEGER :: ij0, ij1 , ikv! temporary integer114 INTEGER :: ii0, ii1 ! temporary integer 115 INTEGER :: ij0, ij1 ! temporary integer 116 116 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_16, zcofw, z1_slpmax ! local scalars 117 117 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - 118 118 REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - - 119 119 REAL(wp) :: zck, zfk, zbw ! - - 120 REAL(wp) :: zdepu, zdepv ! - - 121 REAL(wp), POINTER, DIMENSION(:,: ) :: zslpml_hmlpu, zslpml_hmlpv 120 122 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zww 121 123 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdzr … … 126 128 ! 127 129 CALL wrk_alloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 130 CALL wrk_alloc( jpi,jpj, zslpml_hmlpu, zslpml_hmlpv ) 128 131 129 132 zeps = 1.e-20_wp !== Local constant initialization ==! … … 149 152 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 150 153 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 154 END DO 155 END DO 156 ENDIF 157 IF( ln_zps .AND. ln_isfcav ) THEN ! partial steps correction at the bottom ocean level 158 DO jj = 1, jpjm1 159 DO ji = 1, jpim1 160 IF ( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj) 161 IF ( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj) 151 162 END DO 152 163 END DO … … 171 182 ! =========================== | vslp = d/dj( prd ) / d/dz( prd ) 172 183 ! 184 IF ( ln_isfcav ) THEN 185 DO jj = 2, jpjm1 186 DO ji = fs_2, fs_jpim1 ! vector opt. 187 zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / ( MAX(hmlpt(ji,jj), hmlpt(ji+1,jj ), 5._wp) & 188 & - 0.5_wp * ( risfdep(ji,jj) + risfdep(ji+1,jj ) ) ) 189 zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / ( MAX(hmlpt(ji,jj), hmlpt(ji ,jj+1), 5._wp) & 190 & - 0.5_wp * ( risfdep(ji,jj) + risfdep(ji ,jj+1) ) ) 191 END DO 192 END DO 193 ELSE 194 DO jj = 2, jpjm1 195 DO ji = fs_2, fs_jpim1 ! vector opt. 196 zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji+1,jj ), 5._wp) 197 zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji ,jj+1), 5._wp) 198 END DO 199 END DO 200 END IF 201 173 202 DO jk = 2, jpkm1 !* Slopes at u and v points 174 203 DO jj = 2, jpjm1 … … 186 215 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 187 216 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 188 zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps ) & 189 & + zfi * uslpml(ji,jj) & 190 & * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) ) & 191 & / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk) 192 zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps ) & 193 & + zfj * vslpml(ji,jj) & 194 & * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) ) & 195 & / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk) 217 ! thickness of water column between surface and level k at u/v point 218 zdepu = 0.5_wp * ( ( fsdept (ji,jj,jk) + fsdept (ji+1,jj ,jk) ) & 219 - ( risfdep(ji,jj) + risfdep(ji+1,jj) ) - fse3u(ji,jj,miku(ji,jj)) ) 220 zdepv = 0.5_wp * ( ( fsdept (ji,jj,jk) + fsdept (ji,jj+1,jk) ) & 221 - ( risfdep(ji,jj) + risfdep(ji,jj+1) ) - fse3v(ji,jj,mikv(ji,jj)) ) 222 ! 223 zwz(ji,jj,jk) = ( ( 1._wp - zfi) * zau / ( zbu - zeps ) & 224 & + zfi * zdepu * zslpml_hmlpu(ji,jj) ) * umask(ji,jj,jk) 225 zww(ji,jj,jk) = ( ( 1._wp - zfj) * zav / ( zbv - zeps ) & 226 & + zfj * zdepv * zslpml_hmlpv(ji,jj) ) * vmask(ji,jj,jk) 196 227 !!gm modif to suppress omlmask.... (as in Griffies case) 197 228 ! ! ! jk must be >= ML level for zf=1. otherwise zf=0. … … 265 296 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj) 266 297 zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) & 267 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * tmask (ji,jj,jk)298 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * wmask (ji,jj,jk) 268 299 zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) & 269 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * tmask (ji,jj,jk)300 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * wmask (ji,jj,jk) 270 301 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 271 302 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) … … 274 305 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 275 306 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 276 zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp )277 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * tmask(ji,jj,jk)278 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * tmask(ji,jj,jk)307 zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj) - fsdepw(ji,jj,mikt(ji,jj)), 10._wp ) 308 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * wmask(ji,jj,jk) 309 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * wmask(ji,jj,jk) 279 310 280 311 !!gm modif to suppress omlmask.... (as in Griffies operator) … … 340 371 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 341 372 342 343 373 IF(ln_ctl) THEN 344 374 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) … … 347 377 ! 348 378 CALL wrk_dealloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 379 CALL wrk_dealloc( jpi,jpj, zslpml_hmlpu, zslpml_hmlpv ) 349 380 ! 350 381 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp') … … 486 517 ! 487 518 jk = nmln(ji,jj+jp) + 1 488 IF( jk .GT.mbkt(ji,jj+jp) ) THEN !ML reaches bottom519 IF( jk > mbkt(ji,jj+jp) ) THEN !ML reaches bottom 489 520 ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp 490 521 ELSE … … 699 730 zcj = MAX( vmask(ji,jj-1,ik ) + vmask(ji,jj,ik ) & 700 731 & + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps ) * e2t(ji,jj) 701 zai = ( p_gru(ji-1,jj,ik ) + p_gru(ji,jj,ik) &732 zai = ( p_gru(ji-1,jj,ik ) + p_gru(ji,jj,ik) & 702 733 & + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1 ) ) / zci * tmask(ji,jj,ik) 703 734 zaj = ( p_grv(ji,jj-1,ik ) + p_grv(ji,jj,ik ) & -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
r5836 r6012 44 44 LOGICAL , PUBLIC :: ln_dm2dc !: Daily mean to Diurnal Cycle short wave (qsr) 45 45 LOGICAL , PUBLIC :: ln_rnf !: runoffs / runoff mouths 46 LOGICAL , PUBLIC :: ln_isf !: ice shelf melting 46 47 LOGICAL , PUBLIC :: ln_ssr !: Sea Surface restoring on SST and/or SSS 47 48 LOGICAL , PUBLIC :: ln_apr_dyn !: Atmospheric pressure forcing used on dynamics (ocean & ice) 48 49 INTEGER , PUBLIC :: nn_ice !: flag for ice in the surface boundary condition (=0/1/2/3) 49 INTEGER , PUBLIC :: nn_isf !: flag for isf in the surface boundary condition (=0/1/2/3/4)50 50 INTEGER , PUBLIC :: nn_ice_embd !: flag for levitating/embedding sea-ice in the ocean 51 51 ! !: =0 levitating ice (no mass exchange, concentration/dilution effect) -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
r6006 r6012 35 35 ! public in order to be able to output then 36 36 37 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_tsc_b, risf_tsc 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qisf !: net heat flux from ice shelf37 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_tsc_b, risf_tsc !: before and now T & S isf contents [K.m/s & PSU.m/s] 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qisf !: net heat flux from ice shelf [W/m2] 39 39 REAL(wp), PUBLIC :: rn_hisf_tbl !: thickness of top boundary layer [m] 40 LOGICAL , PUBLIC :: ln_divisf !: flag to correct divergence 41 INTEGER , PUBLIC :: nn_isfblk !: 42 INTEGER , PUBLIC :: nn_gammablk !: 43 LOGICAL , PUBLIC :: ln_conserve !: 44 REAL(wp), PUBLIC :: rn_gammat0 !: temperature exchange coeficient 45 REAL(wp), PUBLIC :: rn_gammas0 !: salinity exchange coeficient 46 REAL(wp), PUBLIC :: rdivisf !: flag to test if fwf apply on divergence 40 INTEGER , PUBLIC :: nn_isf !: flag to choose between explicit/param/specified 41 INTEGER , PUBLIC :: nn_isfblk !: flag to choose the bulk formulation to compute the ice shelf melting 42 INTEGER , PUBLIC :: nn_gammablk !: flag to choose how the exchange coefficient is computed 43 REAL(wp), PUBLIC :: rn_gammat0 !: temperature exchange coeficient [] 44 REAL(wp), PUBLIC :: rn_gammas0 !: salinity exchange coeficient [] 47 45 48 46 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: rzisf_tbl !:depth of calving front (shallowest point) nn_isf ==2/3 49 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: rhisf_tbl, rhisf_tbl_0 !:thickness of tbl 47 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: rhisf_tbl, rhisf_tbl_0 !:thickness of tbl [m] 50 48 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: r1_hisf_tbl !:1/thickness of tbl 51 49 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: ralpha !:proportion of bottom cell influenced by tbl 52 50 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: risfLeff !:effective length (Leff) BG03 nn_isf==2 53 51 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point 54 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: misfkt, misfkb !:Level of ice shelf base55 56 REAL(wp), PUBLIC, SAVE :: rcpi = 2000.0_wp ! phycst ?57 REAL(wp), PUBLIC, SAVE :: kappa = 1.54e-6_wp ! phycst ?58 REAL(wp), PUBLIC, SAVE :: rhoisf = 920.0_wp ! phycst ?59 REAL(wp), PUBLIC, SAVE :: tsurf = -20.0_wp ! phycst ?60 REAL(wp), PUBLIC, SAVE :: lfusisf= 0.334e6_wp ! phycst ?52 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: misfkt, misfkb !:Level of ice shelf base 53 54 REAL(wp), PUBLIC, SAVE :: rcpi = 2000.0_wp ! specific heat of ice shelf [J/kg/K] 55 REAL(wp), PUBLIC, SAVE :: rkappa = 1.54e-6_wp ! heat diffusivity through the ice-shelf [m2/s] 56 REAL(wp), PUBLIC, SAVE :: rhoisf = 920.0_wp ! volumic mass of ice shelf [kg/m3] 57 REAL(wp), PUBLIC, SAVE :: tsurf = -20.0_wp ! air temperature on top of ice shelf [C] 58 REAL(wp), PUBLIC, SAVE :: rlfusisf = 0.334e6_wp ! latent heat of fusion of ice shelf [J/kg] 61 59 62 60 !: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3) 63 CHARACTER(len=100), PUBLIC :: cn_dirisf = './' !: Root directory for location of ssr files 64 TYPE(FLD_N) , PUBLIC :: sn_qisf, sn_fwfisf !: information about the runoff file to be read 65 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_qisf, sf_fwfisf 66 TYPE(FLD_N) , PUBLIC :: sn_rnfisf !: information about the runoff file to be read 67 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf 68 TYPE(FLD_N) , PUBLIC :: sn_depmax_isf, sn_depmin_isf, sn_Leff_isf !: information about the runoff file to be read 61 CHARACTER(len=100), PUBLIC :: cn_dirisf = './' !: Root directory for location of ssr files 62 TYPE(FLD_N) , PUBLIC :: sn_fwfisf !: information about the isf melting file to be read 63 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf 64 TYPE(FLD_N) , PUBLIC :: sn_rnfisf !: information about the isf melting param. file to be read 65 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf 66 TYPE(FLD_N) , PUBLIC :: sn_depmax_isf !: information about the grounding line depth file to be read 67 TYPE(FLD_N) , PUBLIC :: sn_depmin_isf !: information about the calving line depth file to be read 68 TYPE(FLD_N) , PUBLIC :: sn_Leff_isf !: information about the effective length file to be read 69 69 70 70 !! * Substitutions … … 77 77 CONTAINS 78 78 79 79 SUBROUTINE sbc_isf(kt) 80 80 !!--------------------------------------------------------------------- 81 !! *** ROUTINE sbc_isf *** 82 !!--------------------------------------------------------------------- 83 INTEGER, INTENT(in) :: kt ! ocean time step 84 ! 85 INTEGER :: ji, jj, jk, ijkmin, inum, ierror 86 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer 87 REAL(wp) :: rmin 88 REAL(wp) :: zhk 89 REAL(wp) :: zt_frz, zpress 90 CHARACTER(len=256) :: cfisf , cvarzisf, cvarhisf ! name for isf file 91 CHARACTER(LEN=256) :: cnameis ! name of iceshelf file 92 CHARACTER (LEN=32) :: cvarLeff ! variable name for efficient Length scale 93 INTEGER :: ios ! Local integer output status for namelist read 94 !! 95 NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, ln_divisf, ln_conserve, rn_gammat0, rn_gammas0, nn_gammablk, & 96 & sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 81 !! *** ROUTINE sbc_isf *** 82 !! 83 !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf 84 !! melting and freezing 85 !! 86 !! ** Method : 4 parameterizations are available according to nn_isf 87 !! nn_isf = 1 : Realistic ice_shelf formulation 88 !! 2 : Beckmann & Goose parameterization 89 !! 3 : Specified runoff in deptht (Mathiot & al. ) 90 !! 4 : specified fwf and heat flux forcing beneath the ice shelf 91 !!---------------------------------------------------------------------- 92 INTEGER, INTENT( in ) :: kt ! ocean time step 93 ! 94 INTEGER :: ji, jj ! loop index 95 REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep) 97 96 !!--------------------------------------------------------------------- 98 97 ! … … 100 99 IF( kt == nit000 ) THEN ! First call kt=nit000 ! 101 100 ! ! ====================== ! 102 REWIND( numnam_ref ) ! Namelist namsbc_rnf in reference namelist : Runoffs 103 READ ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901) 104 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp ) 105 106 REWIND( numnam_cfg ) ! Namelist namsbc_rnf in configuration namelist : Runoffs 107 READ ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 ) 108 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp ) 109 IF(lwm) WRITE ( numond, namsbc_isf ) 110 111 IF ( lwp ) WRITE(numout,*) 112 IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf' 113 IF ( lwp ) WRITE(numout,*) '~~~~~~~~~' 114 IF ( lwp ) WRITE(numout,*) 'sbcisf :' 115 IF ( lwp ) WRITE(numout,*) '~~~~~~~~' 116 IF ( lwp ) WRITE(numout,*) ' nn_isf = ', nn_isf 117 IF ( lwp ) WRITE(numout,*) ' nn_isfblk = ', nn_isfblk 118 IF ( lwp ) WRITE(numout,*) ' rn_hisf_tbl = ', rn_hisf_tbl 119 IF ( lwp ) WRITE(numout,*) ' ln_divisf = ', ln_divisf 120 IF ( lwp ) WRITE(numout,*) ' nn_gammablk = ', nn_gammablk 121 IF ( lwp ) WRITE(numout,*) ' rn_gammat0 = ', rn_gammat0 122 IF ( lwp ) WRITE(numout,*) ' rn_gammas0 = ', rn_gammas0 123 IF ( lwp ) WRITE(numout,*) ' rn_tfri2 = ', rn_tfri2 124 IF (ln_divisf) THEN ! keep it in the namelist ??? used true anyway as for runoff ? (PM) 125 rdivisf = 1._wp 126 ELSE 127 rdivisf = 0._wp 128 END IF 129 ! 130 ! Allocate public variable 131 IF ( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' ) 132 ! 133 ! initialisation 134 qisf(:,:) = 0._wp ; fwfisf(:,:) = 0._wp 135 risf_tsc(:,:,:) = 0._wp 136 ! 137 ! define isf tbl tickness, top and bottom indice 138 IF (nn_isf == 1) THEN 139 rhisf_tbl(:,:) = rn_hisf_tbl 140 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 141 ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN 142 ALLOCATE( sf_rnfisf(1), STAT=ierror ) 143 ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 144 CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 145 146 !: read effective lenght (BG03) 147 IF (nn_isf == 2) THEN 148 ! Read Data and save some integral values 149 CALL iom_open( sn_Leff_isf%clname, inum ) 150 cvarLeff = 'soLeff' !: variable name for Efficient Length scale 151 CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1) 152 CALL iom_close(inum) 153 ! 154 risfLeff = risfLeff*1000 !: convertion in m 155 END IF 156 157 ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth) 158 CALL iom_open( sn_depmax_isf%clname, inum ) 159 cvarhisf = TRIM(sn_depmax_isf%clvar) 160 CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base 161 CALL iom_close(inum) 162 ! 163 CALL iom_open( sn_depmin_isf%clname, inum ) 164 cvarzisf = TRIM(sn_depmin_isf%clvar) 165 CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base 166 CALL iom_close(inum) 167 ! 168 rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:) !: tickness isf boundary layer 169 170 !! compute first level of the top boundary layer 171 DO ji = 1, jpi 172 DO jj = 1, jpj 173 jk = 2 174 DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdepw(ji,jj,jk) < rzisf_tbl(ji,jj) ) ; jk = jk + 1 ; END DO 175 misfkt(ji,jj) = jk-1 176 END DO 177 END DO 178 179 ELSE IF ( nn_isf == 4 ) THEN 180 ! as in nn_isf == 1 181 rhisf_tbl(:,:) = rn_hisf_tbl 182 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 183 184 ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 185 ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror ) 186 ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 187 ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) ) 188 CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 189 !CALL fld_fill( sf_qisf , (/ sn_qisf /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data' , 'namsbc_isf' ) 190 END IF 191 192 ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 193 rhisf_tbl_0(:,:) = rhisf_tbl(:,:) 194 DO jj = 1,jpj 195 DO ji = 1,jpi 196 ikt = misfkt(ji,jj) 197 ikb = misfkt(ji,jj) 198 ! thickness of boundary layer at least the top level thickness 199 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt)) 200 201 ! determine the deepest level influenced by the boundary layer 202 ! test on tmask useless ????? 203 DO jk = ikt, mbkt(ji,jj) 204 IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 205 END DO 206 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 207 misfkb(ji,jj) = ikb ! last wet level of the tbl 208 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 209 210 zhk = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 211 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer 212 END DO 213 END DO 101 CALL sbc_isf_init 102 ! ! ---------------------------------------- ! 103 ELSE ! Swap of forcing fields ! 104 ! ! ---------------------------------------- ! 105 fwfisf_b (:,: ) = fwfisf (:,: ) ! Swap the ocean forcing fields except at nit000 106 risf_tsc_b(:,:,:) = risf_tsc(:,:,:) ! where before fields are set at the end of the routine 214 107 ! 215 108 END IF 216 109 217 ! ! ---------------------------------------- !218 IF( kt /= nit000 ) THEN ! Swap of forcing fields !219 ! ! ---------------------------------------- !220 fwfisf_b (:,: ) = fwfisf (:,: ) ! Swap the ocean forcing fields except at nit000221 risf_tsc_b(:,:,:) = risf_tsc(:,:,:) ! where before fields are set at the end of the routine222 !223 ENDIF224 225 110 IF( MOD( kt-1, nn_fsbc) == 0 ) THEN 226 227 228 ! compute salf and heat flux 229 IF (nn_isf == 1) THEN 230 ! realistic ice shelf formulation 111 ! allocation 112 CALL wrk_alloc( jpi,jpj, zt_frz, zdep ) 113 114 ! compute salt and heat flux 115 SELECT CASE ( nn_isf ) 116 CASE ( 1 ) ! realistic ice shelf formulation 231 117 ! compute T/S/U/V for the top boundary layer 232 118 CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') … … 242 128 CALL sbc_isf_cav (kt) 243 129 244 ELSE IF (nn_isf == 2) THEN 245 ! Beckmann and Goosse parametrisation 130 CASE ( 2 ) ! Beckmann and Goosse parametrisation 246 131 stbl(:,:) = soce 247 132 CALL sbc_isf_bg03(kt) 248 133 249 ELSE IF (nn_isf == 3) THEN 250 ! specified runoff in depth (Mathiot et al., XXXX in preparation) 134 CASE ( 3 ) ! specified runoff in depth (Mathiot et al., XXXX in preparation) 251 135 CALL fld_read ( kt, nn_fsbc, sf_rnfisf ) 252 fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1) ! f resh waterflux from the isf (fwfisf <0 mean melting)253 qisf(:,:) = fwfisf(:,:) * lfusisf ! heatflux136 fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1) ! fwf flux from the isf (fwfisf <0 mean melting) 137 qisf(:,:) = fwfisf(:,:) * rlfusisf ! heat flux 254 138 stbl(:,:) = soce 255 139 256 ELSE IF (nn_isf == 4) THEN 257 ! specified fwf and heat flux forcing beneath the ice shelf 140 CASE ( 4 ) ! specified fwf and heat flux forcing beneath the ice shelf 258 141 CALL fld_read ( kt, nn_fsbc, sf_fwfisf ) 259 !CALL fld_read ( kt, nn_fsbc, sf_qisf ) 260 fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1) ! fwf 261 qisf(:,:) = fwfisf(:,:) * lfusisf ! heat flux 262 !qisf(:,:) = sf_qisf(1)%fnow(:,:,1) ! heat flux 142 fwfisf(:,:) = - sf_fwfisf(1)%fnow(:,:,1) ! fwf flux from the isf (fwfisf <0 mean melting) 143 qisf(:,:) = fwfisf(:,:) * rlfusisf ! heat flux 263 144 stbl(:,:) = soce 264 145 265 END IF 146 END SELECT 147 266 148 ! compute tsc due to isf 267 ! WARNING water add at temp = 0C, correction term is added, maybe better here but need a 3D variable). 268 ! zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 269 zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress ) 270 risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - rdivisf * fwfisf(:,:) * zt_frz * r1_rau0 ! 149 ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU. 150 ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rau0). 151 ! compute freezing point beneath ice shelf (or top cell if nn_isf = 3) 152 DO jj = 1,jpj 153 DO ji = 1,jpi 154 zdep(ji,jj)=fsdepw_n(ji,jj,misfkt(ji,jj)) 155 END DO 156 END DO 157 CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) ) 271 158 272 ! salt effect already take into account in vertical advection 273 risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0 274 275 ! output 276 IF( iom_use('qisf' ) ) CALL iom_put('qisf' , qisf) 277 IF( iom_use('fwfisf') ) CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce ) 278 279 ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now 280 fwfisf(:,:) = rdivisf * fwfisf(:,:) 281 159 risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 ! 160 risf_tsc(:,:,jp_sal) = 0.0_wp 161 282 162 ! lbclnk 283 163 CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.) 284 164 CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.) 285 CALL lbc_lnk(fwfisf(:,:) ,'T',1.)286 CALL lbc_lnk(qisf(:,:) ,'T',1.)287 288 IF( kt == nit000 ) THEN 165 CALL lbc_lnk(fwfisf(:,:) ,'T',1.) 166 CALL lbc_lnk(qisf(:,:) ,'T',1.) 167 168 IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 ! 289 169 IF( ln_rstart .AND. & ! Restart: read in restart file 290 170 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN … … 297 177 risf_tsc_b(:,:,:)= risf_tsc(:,:,:) 298 178 END IF 299 END IF179 END IF 300 180 ! 181 ! output 182 CALL iom_put('qisf' , qisf) 183 CALL iom_put('fwfisf', fwfisf) 184 185 ! deallocation 186 CALL wrk_dealloc( jpi,jpj, zt_frz, zdep ) 301 187 END IF 302 188 ! … … 317 203 & STAT= sbc_isf_alloc ) 318 204 ! 319 IF( lk_mpp 205 IF( lk_mpp ) CALL mpp_sum ( sbc_isf_alloc ) 320 206 IF( sbc_isf_alloc /= 0 ) CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.') 321 207 ! 322 END IF208 END IF 323 209 END FUNCTION 324 210 325 326 SUBROUTINE sbc_isf_bg03(kt) 327 !!========================================================================== 328 !! *** SUBROUTINE sbcisf_bg03 *** 329 !! add net heat and fresh water flux from ice shelf melting 330 !! into the adjacent ocean using the parameterisation by 331 !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 332 !! interaction for climate models", Ocean Modelling 5(2003) 157-170. 333 !! (hereafter BG) 334 !!========================================================================== 335 !!---------------------------------------------------------------------- 336 !! sbc_isf_bg03 : routine called from sbcmod 337 !!---------------------------------------------------------------------- 338 !! 339 !! ** Purpose : Add heat and fresh water fluxes due to ice shelf melting 340 !! ** Reference : Beckmann et Goosse, 2003, Ocean Modelling 341 !! 211 SUBROUTINE sbc_isf_init 212 !!--------------------------------------------------------------------- 213 !! *** ROUTINE sbc_isf_init *** 214 !! 215 !! ** Purpose : Initialisation of variables for iceshelf fluxes formulation 216 !! 217 !! ** Method : 4 parameterizations are available according to nn_isf 218 !! nn_isf = 1 : Realistic ice_shelf formulation 219 !! 2 : Beckmann & Goose parameterization 220 !! 3 : Specified runoff in deptht (Mathiot & al. ) 221 !! 4 : specified fwf and heat flux forcing beneath the ice shelf 222 !!---------------------------------------------------------------------- 223 INTEGER :: ji, jj, jk ! loop index 224 INTEGER :: ik ! current level index 225 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer 226 INTEGER :: inum, ierror 227 INTEGER :: ios ! Local integer output status for namelist read 228 REAL(wp) :: zhk 229 CHARACTER(len=256) :: cvarzisf, cvarhisf ! name for isf file 230 CHARACTER(LEN=32 ) :: cvarLeff ! variable name for efficient Length scale 231 !!---------------------------------------------------------------------- 232 NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf, & 233 & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 234 !!---------------------------------------------------------------------- 235 236 REWIND( numnam_ref ) ! Namelist namsbc_rnf in reference namelist : Runoffs 237 READ ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901) 238 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp ) 239 240 REWIND( numnam_cfg ) ! Namelist namsbc_rnf in configuration namelist : Runoffs 241 READ ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 ) 242 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp ) 243 IF(lwm) WRITE ( numond, namsbc_isf ) 244 245 IF ( lwp ) WRITE(numout,*) 246 IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf' 247 IF ( lwp ) WRITE(numout,*) '~~~~~~~~~' 248 IF ( lwp ) WRITE(numout,*) 'sbcisf :' 249 IF ( lwp ) WRITE(numout,*) '~~~~~~~~' 250 IF ( lwp ) WRITE(numout,*) ' nn_isf = ', nn_isf 251 IF ( lwp ) WRITE(numout,*) ' nn_isfblk = ', nn_isfblk 252 IF ( lwp ) WRITE(numout,*) ' rn_hisf_tbl = ', rn_hisf_tbl 253 IF ( lwp ) WRITE(numout,*) ' nn_gammablk = ', nn_gammablk 254 IF ( lwp ) WRITE(numout,*) ' rn_gammat0 = ', rn_gammat0 255 IF ( lwp ) WRITE(numout,*) ' rn_gammas0 = ', rn_gammas0 256 IF ( lwp ) WRITE(numout,*) ' rn_tfri2 = ', rn_tfri2 257 ! 258 ! Allocate public variable 259 IF ( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' ) 260 ! 261 ! initialisation 262 qisf(:,:) = 0._wp ; fwfisf (:,:) = 0._wp 263 risf_tsc(:,:,:) = 0._wp ; fwfisf_b(:,:) = 0._wp 264 ! 265 ! define isf tbl tickness, top and bottom indice 266 SELECT CASE ( nn_isf ) 267 CASE ( 1 ) 268 rhisf_tbl(:,:) = rn_hisf_tbl 269 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 270 271 CASE ( 2 , 3 ) 272 ALLOCATE( sf_rnfisf(1), STAT=ierror ) 273 ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 274 CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 275 276 ! read effective lenght (BG03) 277 IF (nn_isf == 2) THEN 278 CALL iom_open( sn_Leff_isf%clname, inum ) 279 cvarLeff = TRIM(sn_Leff_isf%clvar) 280 CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1) 281 CALL iom_close(inum) 282 ! 283 risfLeff = risfLeff*1000.0_wp !: convertion in m 284 END IF 285 286 ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth) 287 CALL iom_open( sn_depmax_isf%clname, inum ) 288 cvarhisf = TRIM(sn_depmax_isf%clvar) 289 CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base 290 CALL iom_close(inum) 291 ! 292 CALL iom_open( sn_depmin_isf%clname, inum ) 293 cvarzisf = TRIM(sn_depmin_isf%clvar) 294 CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base 295 CALL iom_close(inum) 296 ! 297 rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:) !: tickness isf boundary layer 298 299 !! compute first level of the top boundary layer 300 DO ji = 1, jpi 301 DO jj = 1, jpj 302 ik = 2 303 DO WHILE ( ik <= mbkt(ji,jj) .AND. fsdepw(ji,jj,ik) < rzisf_tbl(ji,jj) ) ; ik = ik + 1 ; END DO 304 misfkt(ji,jj) = ik-1 305 END DO 306 END DO 307 308 CASE ( 4 ) 309 ! as in nn_isf == 1 310 rhisf_tbl(:,:) = rn_hisf_tbl 311 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 312 313 ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 314 ALLOCATE( sf_fwfisf(1), STAT=ierror ) 315 ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 316 CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 317 318 END SELECT 319 320 rhisf_tbl_0(:,:) = rhisf_tbl(:,:) 321 322 ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 323 DO jj = 1,jpj 324 DO ji = 1,jpi 325 ikt = misfkt(ji,jj) 326 ikb = misfkt(ji,jj) 327 ! thickness of boundary layer at least the top level thickness 328 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt)) 329 330 ! determine the deepest level influenced by the boundary layer 331 DO jk = ikt+1, mbkt(ji,jj) 332 IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 333 END DO 334 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 335 misfkb(ji,jj) = ikb ! last wet level of the tbl 336 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 337 338 zhk = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 339 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer 340 END DO 341 END DO 342 343 END SUBROUTINE sbc_isf_init 344 345 SUBROUTINE sbc_isf_bg03(kt) 346 !!--------------------------------------------------------------------- 347 !! *** ROUTINE sbc_isf_bg03 *** 348 !! 349 !! ** Purpose : add net heat and fresh water flux from ice shelf melting 350 !! into the adjacent ocean 351 !! 352 !! ** Method : See reference 353 !! 354 !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 355 !! interaction for climate models", Ocean Modelling 5(2003) 157-170. 356 !! (hereafter BG) 342 357 !! History : 343 !! !06-02 (C. Wang) Original code358 !! 06-02 (C. Wang) Original code 344 359 !!---------------------------------------------------------------------- 345 360 INTEGER, INTENT ( in ) :: kt 346 361 ! 347 INTEGER :: ji, jj, jk, jish !temporary integer 348 INTEGER :: ijkmin 349 INTEGER :: ii, ij, ik 350 INTEGER :: inum 351 352 REAL(wp) :: zt_sum ! sum of the temperature between 200m and 600m 353 REAL(wp) :: zt_ave ! averaged temperature between 200m and 600m 354 REAL(wp) :: zt_frz ! freezing point temperature at depth z 355 REAL(wp) :: zpress ! pressure to compute the freezing point in depth 356 357 !!---------------------------------------------------------------------- 358 IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 359 ! 360 DO ji = 1, jpi 361 DO jj = 1, jpj 362 ik = misfkt(ji,jj) 363 !! Initialize arrays to 0 (each step) 364 zt_sum = 0.e0_wp 365 IF ( ik .GT. 1 ) THEN 366 ! 3. -----------the average temperature between 200m and 600m --------------------- 367 DO jk = misfkt(ji,jj),misfkb(ji,jj) 368 ! Calculate freezing temperature 369 CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, fsdept(ji,jj,ik)) 370 zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik) ! sum temp 371 ENDDO 372 zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 373 374 ! 4. ------------Net heat flux and fresh water flux due to the ice shelf 375 ! For those corresponding to zonal boundary 376 qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave & 377 & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,ik) 362 INTEGER :: ji, jj, jk ! dummy loop index 363 INTEGER :: ik ! current level 364 REAL(wp) :: zt_sum ! sum of the temperature between 200m and 600m 365 REAL(wp) :: zt_ave ! averaged temperature between 200m and 600m 366 REAL(wp) :: zt_frz ! freezing point temperature at depth z 367 REAL(wp) :: zpress ! pressure to compute the freezing point in depth 368 !!---------------------------------------------------------------------- 369 370 IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 371 ! 372 DO ji = 1, jpi 373 DO jj = 1, jpj 374 ik = misfkt(ji,jj) 375 !! Initialize arrays to 0 (each step) 376 zt_sum = 0.e0_wp 377 IF ( ik > 1 ) THEN 378 ! 1. -----------the average temperature between 200m and 600m --------------------- 379 DO jk = misfkt(ji,jj),misfkb(ji,jj) 380 ! freezing point temperature at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!) 381 ! after verif with UNESCO, wrong sign in BG eq. 2 382 ! Calculate freezing temperature 383 CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 384 zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) ! sum temp 385 END DO 386 zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 387 ! 2. ------------Net heat flux and fresh water flux due to the ice shelf 388 ! For those corresponding to zonal boundary 389 qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave & 390 & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk) 378 391 379 fwfisf(ji,jj) = qisf(ji,jj) /lfusisf !fresh water flux kg/(m2s)380 fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )381 !add to salinity trend382 ELSE383 qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp384 END IF385 END DO386 END DO387 !388 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_bg03')392 fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf !fresh water flux kg/(m2s) 393 fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 394 !add to salinity trend 395 ELSE 396 qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 397 END IF 398 END DO 399 END DO 400 ! 401 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_bg03') 389 402 ! 390 403 END SUBROUTINE sbc_isf_bg03 391 404 392 393 SUBROUTINE sbc_isf_cav( kt ) 405 SUBROUTINE sbc_isf_cav( kt ) 394 406 !!--------------------------------------------------------------------- 395 407 !! *** ROUTINE sbc_isf_cav *** … … 406 418 INTEGER, INTENT(in) :: kt ! ocean time step 407 419 ! 408 LOGICAL :: ln_isomip = .true. 409 REAL(wp), DIMENSION(:,:), POINTER :: zfrz,zpress,zti 410 REAL(wp), DIMENSION(:,:), POINTER :: zgammat2d, zgammas2d 411 !REAL(wp), DIMENSION(:,:), POINTER :: zqisf, zfwfisf 420 INTEGER :: ji, jj ! dummy loop indices 421 INTEGER :: nit 412 422 REAL(wp) :: zlamb1, zlamb2, zlamb3 413 423 REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7 414 424 REAL(wp) :: zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac 415 REAL(wp) :: zfwflx, zhtflx, zhtflx_b 416 REAL(wp) :: zgammat, zgammas 417 REAL(wp) :: zeps = 1.e-20_wp !== Local constant initialization ==! 418 INTEGER :: ji, jj ! dummy loop indices 419 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers 420 INTEGER :: ierror ! return error code 421 LOGICAL :: lit=.TRUE. 422 INTEGER :: nit 425 REAL(wp) :: zeps = 1.e-20_wp 426 REAL(wp) :: zerr 427 REAL(wp), DIMENSION(:,:), POINTER :: zfrz 428 REAL(wp), DIMENSION(:,:), POINTER :: zgammat, zgammas 429 REAL(wp), DIMENSION(:,:), POINTER :: zfwflx, zhtflx, zhtflx_b 430 LOGICAL :: lit 423 431 !!--------------------------------------------------------------------- 424 ! 425 ! coeficient for linearisation of tfreez426 zlamb1 =-0.0575427 zlamb2 =0.0901428 zlamb3 =-7.61e-04432 ! coeficient for linearisation of potential tfreez 433 ! Crude approximation for pressure (but commonly used) 434 zlamb1 =-0.0573_wp 435 zlamb2 = 0.0832_wp 436 zlamb3 =-7.53e-08_wp * grav * rau0 429 437 IF( nn_timing == 1 ) CALL timing_start('sbc_isf_cav') 430 438 ! 431 CALL wrk_alloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d ) 432 433 zcfac=0.0_wp 434 IF (ln_conserve) zcfac=1.0_wp 435 zpress(:,:)=0.0_wp 436 zgammat2d(:,:)=0.0_wp 437 zgammas2d(:,:)=0.0_wp 438 ! 439 ! 440 DO jj = 1, jpj 441 DO ji = 1, jpi 442 ! Crude approximation for pressure (but commonly used) 443 ! 1e-04 to convert from Pa to dBar 444 zpress(ji,jj)=grav*rau0*fsdepw(ji,jj,mikt(ji,jj))*1.e-04 445 ! 446 END DO 439 CALL wrk_alloc( jpi,jpj, zfrz , zgammat, zgammas ) 440 CALL wrk_alloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b ) 441 442 ! initialisation 443 zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0 444 zhtflx (:,:) = 0.0_wp ; zhtflx_b(:,:) = 0.0_wp 445 zfwflx (:,:) = 0.0_wp 446 447 ! compute ice shelf melting 448 nit = 1 ; lit = .TRUE. 449 DO WHILE ( lit ) ! maybe just a constant number of iteration as in blk_core is fine 450 SELECT CASE ( nn_isfblk ) 451 CASE ( 1 ) ! ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006) 452 ! Calculate freezing temperature 453 CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) ) 454 455 ! compute gammat every where (2d) 456 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 457 458 ! compute upward heat flux zhtflx and upward water flux zwflx 459 DO jj = 1, jpj 460 DO ji = 1, jpi 461 zhtflx(ji,jj) = zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj)) 462 zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf 463 END DO 464 END DO 465 466 ! Compute heat flux and upward fresh water flux 467 qisf (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 468 fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 469 470 CASE ( 2 ) ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 471 ! compute gammat every where (2d) 472 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 473 474 ! compute upward heat flux zhtflx and upward water flux zwflx 475 ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015) 476 DO jj = 1, jpj 477 DO ji = 1, jpi 478 ! compute coeficient to solve the 2nd order equation 479 zeps1 = rcp*rau0*zgammat(ji,jj) 480 zeps2 = rlfusisf*rau0*zgammas(ji,jj) 481 zeps3 = rhoisf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps) 482 zeps4 = zlamb2+zlamb3*risfdep(ji,jj) 483 zeps6 = zeps4-ttbl(ji,jj) 484 zeps7 = zeps4-tsurf 485 zaqe = zlamb1 * (zeps1 + zeps3) 486 zaqer = 0.5_wp/MIN(zaqe,-zeps) 487 zbqe = zeps1*zeps6+zeps3*zeps7-zeps2 488 zcqe = zeps2*stbl(ji,jj) 489 zdis = zbqe*zbqe-4.0_wp*zaqe*zcqe 490 491 ! Presumably zdis can never be negative because gammas is very small compared to gammat 492 ! compute s freeze 493 zsfrz=(-zbqe-SQRT(zdis))*zaqer 494 IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer 495 496 ! compute t freeze (eq. 22) 497 zfrz(ji,jj)=zeps4+zlamb1*zsfrz 498 499 ! zfwflx is upward water flux 500 ! zhtflx is upward heat flux (out of ocean) 501 ! compute the upward water and heat flux (eq. 28 and eq. 29) 502 zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) 503 zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) ) 504 END DO 505 END DO 506 507 ! compute heat and water flux 508 qisf (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 509 fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 510 511 END SELECT 512 513 ! define if we need to iterate (nn_gammablk 0/1 do not need iteration) 514 IF ( nn_gammablk < 2 ) THEN ; lit = .FALSE. 515 ELSE 516 ! check total number of iteration 517 IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 518 ELSE ; nit = nit + 1 519 END IF 520 521 ! compute error between 2 iterations 522 ! if needed save gammat and compute zhtflx_b for next iteration 523 zerr = MAXVAL(ABS(zhtflx-zhtflx_b)) 524 IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE. 525 ELSE ; zhtflx_b(:,:) = zhtflx(:,:) 526 END IF 527 END IF 447 528 END DO 448 449 ! Calculate in-situ temperature (ref to surface) 450 zti(:,:)=tinsitu( ttbl, stbl, zpress ) 451 ! Calculate freezing temperature 452 CALL eos_fzp( sss_m(:,:), zfrz(:,:), zpress ) 453 454 455 zhtflx=0._wp ; zfwflx=0._wp 456 IF (nn_isfblk == 1) THEN 457 DO jj = 1, jpj 458 DO ji = 1, jpi 459 IF (mikt(ji,jj) > 1 ) THEN 460 nit = 1; lit = .TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp 461 DO WHILE ( lit ) 462 ! compute gamma 463 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit) 464 ! zhtflx is upward heat flux (out of ocean) 465 zhtflx = zgammat*rcp*rau0*(zti(ji,jj)-zfrz(ji,jj)) 466 ! zwflx is upward water flux 467 zfwflx = - zhtflx/lfusisf 468 ! test convergence and compute gammat 469 IF ( (zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE. 470 471 nit = nit + 1 472 IF (nit .GE. 100) CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 473 474 ! save gammat and compute zhtflx_b 475 zgammat2d(ji,jj)=zgammat 476 zhtflx_b = zhtflx 477 END DO 478 479 qisf(ji,jj) = - zhtflx 480 ! For genuine ISOMIP protocol this should probably be something like 481 fwfisf(ji,jj) = zfwflx * ( soce / MAX(stbl(ji,jj),zeps)) 482 ELSE 483 fwfisf(ji,jj) = 0._wp 484 qisf(ji,jj) = 0._wp 485 END IF 486 ! 487 END DO 488 END DO 489 490 ELSE IF (nn_isfblk == 2 ) THEN 491 492 ! More complicated 3 equation thermodynamics as in MITgcm 493 DO jj = 2, jpj 494 DO ji = 2, jpi 495 IF (mikt(ji,jj) > 1 ) THEN 496 nit=1; lit=.TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp; zhtflx=0._wp 497 DO WHILE ( lit ) 498 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit) 499 500 zeps1=rcp*rau0*zgammat 501 zeps2=lfusisf*rau0*zgammas 502 zeps3=rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps) 503 zeps4=zlamb2+zlamb3*risfdep(ji,jj) 504 zeps6=zeps4-zti(ji,jj) 505 zeps7=zeps4-tsurf 506 zaqe=zlamb1 * (zeps1 + zeps3) 507 zaqer=0.5/MIN(zaqe,-zeps) 508 zbqe=zeps1*zeps6+zeps3*zeps7-zeps2 509 zcqe=zeps2*stbl(ji,jj) 510 zdis=zbqe*zbqe-4.0*zaqe*zcqe 511 ! Presumably zdis can never be negative because gammas is very small compared to gammat 512 zsfrz=(-zbqe-SQRT(zdis))*zaqer 513 IF (zsfrz .lt. 0.0) zsfrz=(-zbqe+SQRT(zdis))*zaqer 514 zfrz(ji,jj)=zeps4+zlamb1*zsfrz 515 516 ! zfwflx is upward water flux (MAX usefull if kappa = 0 517 zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) ) 518 ! zhtflx is upward heat flux (out of ocean) 519 ! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value 520 zhtflx = ( zgammat*rau0 - zcfac*zfwflx ) * rcp * (zti(ji,jj) - zfrz(ji,jj) ) 521 ! zwflx is upward water flux 522 ! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz 523 zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj) 524 ! test convergence and compute gammat 525 IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE. 526 527 nit = nit + 1 528 IF (nit .GE. 51) THEN 529 WRITE(numout,*) "sbcisf : too many iteration ... ", & 530 & zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea 531 CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 532 END IF 533 ! save gammat and compute zhtflx_b 534 zgammat2d(ji,jj)=zgammat 535 zgammas2d(ji,jj)=zgammas 536 zhtflx_b = zhtflx 537 538 END DO 539 ! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value 540 qisf(ji,jj) = - zhtflx 541 ! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz 542 fwfisf(ji,jj) = zfwflx 543 ELSE 544 fwfisf(ji,jj) = 0._wp 545 qisf(ji,jj) = 0._wp 546 ENDIF 547 ! 548 END DO 549 END DO 550 ENDIF 551 ! lbclnk 552 CALL lbc_lnk(zgammas2d(:,:),'T',1.) 553 CALL lbc_lnk(zgammat2d(:,:),'T',1.) 554 ! output 555 CALL iom_put('isfgammat', zgammat2d) 556 CALL iom_put('isfgammas', zgammas2d) 557 ! 558 CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d ) 529 ! 530 CALL iom_put('isfgammat', zgammat) 531 CALL iom_put('isfgammas', zgammas) 532 ! 533 CALL wrk_dealloc( jpi,jpj, zfrz , zgammat, zgammas ) 534 CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b ) 559 535 ! 560 536 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_cav') … … 562 538 END SUBROUTINE sbc_isf_cav 563 539 564 565 SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit ) 540 SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf ) 566 541 !!---------------------------------------------------------------------- 567 542 !! ** Purpose : compute the coefficient echange for heat flux … … 572 547 !! Jenkins et al., 2010, JPO, p2298-2312 573 548 !!--------------------------------------------------------------------- 574 REAL(wp), INTENT(inout) :: gt, gs, zqhisf, zqwisf575 INTEGER , INTENT(in) :: ji,jj576 LOGICAL , INTENT(inout) :: lit577 578 INTEGER :: ikt! loop index579 REAL(wp) :: zut, zvt,zustar ! U, V at T point and friction velocity549 REAL(wp), DIMENSION(:,:), INTENT(out) :: pgt, pgs 550 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pqhisf, pqwisf 551 ! 552 INTEGER :: ikt 553 INTEGER :: ji, jj ! loop index 554 REAL(wp), DIMENSION(:,:), POINTER :: zustar ! U, V at T point and friction velocity 580 555 REAL(wp) :: zdku, zdkv ! U, V shear 581 556 REAL(wp) :: zPr, zSc, zRc ! Prandtl, Scmidth and Richardson number … … 587 562 REAL(wp) :: zcoef ! temporary coef 588 563 REAL(wp) :: zdep 589 REAL(wp), PARAMETER :: zxsiN = 0.052 ! dimensionless constant 590 REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number 591 REAL(wp), PARAMETER :: znu = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1) 592 REAL(wp) :: rcs = 1.0e-3_wp ! conversion: mm/s ==> m/s 564 REAL(wp) :: zeps = 1.0e-20_wp 565 REAL(wp), PARAMETER :: zxsiN = 0.052_wp ! dimensionless constant 566 REAL(wp), PARAMETER :: znu = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1) 593 567 REAL(wp), DIMENSION(2) :: zts, zab 594 568 !!--------------------------------------------------------------------- 595 ! 596 IF( nn_gammablk == 0 ) THEN 597 !! gamma is constant (specified in namelist) 598 gt = rn_gammat0 599 gs = rn_gammas0 600 lit = .FALSE. 601 ELSE IF ( nn_gammablk == 1 ) THEN 602 !! gamma is assume to be proportional to u* 603 !! WARNING in case of Losh 2008 tbl parametrization, 604 !! you have to used the mean value of u in the boundary layer) 605 !! not yet coded 606 !! Jenkins et al., 2010, JPO, p2298-2312 607 ikt = mikt(ji,jj) 608 !! Compute U and V at T points 609 ! zut = 0.5 * ( utbl(ji-1,jj ) + utbl(ji,jj) ) 610 ! zvt = 0.5 * ( vtbl(ji ,jj-1) + vtbl(ji,jj) ) 611 zut = utbl(ji,jj) 612 zvt = vtbl(ji,jj) 613 614 !! compute ustar 615 zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) ) 616 !! Compute mean value over the TBL 617 618 !! Compute gammats 619 gt = zustar * rn_gammat0 620 gs = zustar * rn_gammas0 621 lit = .FALSE. 622 ELSE IF ( nn_gammablk == 2 ) THEN 623 !! gamma depends of stability of boundary layer 624 !! WARNING in case of Losh 2008 tbl parametrization, 625 !! you have to used the mean value of u in the boundary layer) 626 !! not yet coded 627 !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 628 !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 569 CALL wrk_alloc( jpi,jpj, zustar ) 570 ! 571 SELECT CASE ( nn_gammablk ) 572 CASE ( 0 ) ! gamma is constant (specified in namelist) 573 !! ISOMIP formulation (Hunter et al, 2006) 574 pgt(:,:) = rn_gammat0 575 pgs(:,:) = rn_gammas0 576 577 CASE ( 1 ) ! gamma is assume to be proportional to u* 578 !! Jenkins et al., 2010, JPO, p2298-2312 579 !! Adopted by Asay-Davis et al. (2015) 580 581 !! compute ustar (eq. 24) 582 zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 583 584 !! Compute gammats 585 pgt(:,:) = zustar(:,:) * rn_gammat0 586 pgs(:,:) = zustar(:,:) * rn_gammas0 587 588 CASE ( 2 ) ! gamma depends of stability of boundary layer 589 !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 590 !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 591 !! compute ustar 592 zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 593 594 !! compute Pr and Sc number (can be improved) 595 zPr = 13.8_wp 596 zSc = 2432.0_wp 597 598 !! compute gamma mole 599 zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp 600 zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp 601 602 !! compute gamma 603 DO ji=2,jpi 604 DO jj=2,jpj 629 605 ikt = mikt(ji,jj) 630 606 631 !! Compute U and V at T points 632 zut = 0.5 * ( utbl(ji-1,jj ) + utbl(ji,jj) ) 633 zvt = 0.5 * ( vtbl(ji ,jj-1) + vtbl(ji,jj) ) 634 635 !! compute ustar 636 zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) ) 637 IF (zustar == 0._wp) THEN ! only for kt = 1 I think 638 gt = rn_gammat0 639 gs = rn_gammas0 607 IF (zustar(ji,jj) == 0._wp) THEN ! only for kt = 1 I think 608 pgt = rn_gammat0 609 pgs = rn_gammas0 640 610 ELSE 641 !! compute Rc number (as done in zdfric.F90) 642 zcoef = 0.5 / fse3w(ji,jj,ikt) 643 ! ! shear of horizontal velocity 644 zdku = zcoef * ( un(ji-1,jj ,ikt ) + un(ji,jj,ikt ) & 645 & -un(ji-1,jj ,ikt+1) - un(ji,jj,ikt+1) ) 646 zdkv = zcoef * ( vn(ji ,jj-1,ikt ) + vn(ji,jj,ikt ) & 647 & -vn(ji ,jj-1,ikt+1) - vn(ji,jj,ikt+1) ) 648 ! ! richardson number (minimum value set to zero) 649 zRc = rn2(ji,jj,ikt+1) / ( zdku*zdku + zdkv*zdkv + 1.e-20 ) 650 651 !! compute Pr and Sc number (can be improved) 652 zPr = 13.8 653 zSc = 2432.0 654 655 !! compute gamma mole 656 zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0 657 zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0 658 659 !! compute bouyancy 660 zts(jp_tem) = ttbl(ji,jj) 661 zts(jp_sal) = stbl(ji,jj) 662 zdep = fsdepw(ji,jj,ikt) 663 ! 664 CALL eos_rab( zts, zdep, zab ) 611 !! compute Rc number (as done in zdfric.F90) 612 zcoef = 0.5_wp / fse3w(ji,jj,ikt) 613 ! ! shear of horizontal velocity 614 zdku = zcoef * ( un(ji-1,jj ,ikt ) + un(ji,jj,ikt ) & 615 & -un(ji-1,jj ,ikt+1) - un(ji,jj,ikt+1) ) 616 zdkv = zcoef * ( vn(ji ,jj-1,ikt ) + vn(ji,jj,ikt ) & 617 & -vn(ji ,jj-1,ikt+1) - vn(ji,jj,ikt+1) ) 618 ! ! richardson number (minimum value set to zero) 619 zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps ) 620 621 !! compute bouyancy 622 zts(jp_tem) = ttbl(ji,jj) 623 zts(jp_sal) = stbl(ji,jj) 624 zdep = fsdepw(ji,jj,ikt) 665 625 ! 666 !! compute length scale 667 zbuofdep = grav * ( zab(jp_tem) * zqhisf - zab(jp_sal) * zqwisf ) !!!!!!!!!!!!!!!!!!!!!!!!!!!! 668 669 !! compute Monin Obukov Length 670 ! Maximum boundary layer depth 671 zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) -0.001 672 ! Compute Monin obukhov length scale at the surface and Ekman depth: 673 zmob = zustar ** 3 / (vkarmn * (zbuofdep + epsln)) 674 zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 675 676 !! compute eta* (stability parameter) 677 zetastar = 1 / ( SQRT(1 + MAX(zxsiN * zustar / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0))) 678 679 !! compute the sublayer thickness 680 zhnu = 5 * znu / zustar 681 !! compute gamma turb 682 zgturb = 1/vkarmn * LOG(zustar * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) & 683 & + 1 / ( 2 * zxsiN * zetastar ) - 1/vkarmn 684 685 !! compute gammats 686 gt = zustar / (zgturb + zgmolet) 687 gs = zustar / (zgturb + zgmoles) 626 CALL eos_rab( zts, zdep, zab ) 627 ! 628 !! compute length scale 629 zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) ) !!!!!!!!!!!!!!!!!!!!!!!!!!!! 630 631 !! compute Monin Obukov Length 632 ! Maximum boundary layer depth 633 zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) - 0.001_wp 634 ! Compute Monin obukhov length scale at the surface and Ekman depth: 635 zmob = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps)) 636 zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 637 638 !! compute eta* (stability parameter) 639 zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0_wp))) 640 641 !! compute the sublayer thickness 642 zhnu = 5 * znu / zustar(ji,jj) 643 644 !! compute gamma turb 645 zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) & 646 & + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn 647 648 !! compute gammats 649 pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 650 pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 688 651 END IF 689 END IF 690 ! 691 END SUBROUTINE 692 693 694 SUBROUTINE sbc_isf_tbl( varin, varout, cptin ) 652 END DO 653 END DO 654 CALL lbc_lnk(pgt(:,:),'T',1.) 655 CALL lbc_lnk(pgs(:,:),'T',1.) 656 END SELECT 657 CALL wrk_dealloc( jpi,jpj, zustar ) 658 ! 659 END SUBROUTINE sbc_isf_gammats 660 661 SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin ) 695 662 !!---------------------------------------------------------------------- 696 663 !! *** SUBROUTINE sbc_isf_tbl *** 697 664 !! 698 !! ** Purpose : compute mean T/S/U/V in the boundary layer 699 !! 700 !!---------------------------------------------------------------------- 701 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin 702 REAL(wp), DIMENSION(:,:) , INTENT(out):: varout 665 !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point 666 !! 667 !!---------------------------------------------------------------------- 668 REAL(wp), DIMENSION(:,:,:), INTENT( in ) :: pvarin 669 REAL(wp), DIMENSION(:,:) , INTENT( out ) :: pvarout 670 CHARACTER(len=1), INTENT( in ) :: cd_ptin ! point of variable in/out 671 ! 672 REAL(wp) :: ze3, zhk 673 REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl ! thickness of the tbl 674 675 INTEGER :: ji, jj, jk ! loop index 676 INTEGER :: ikt, ikb ! top and bottom index of the tbl 677 !!---------------------------------------------------------------------- 678 ! allocation 679 CALL wrk_alloc( jpi,jpj, zhisf_tbl) 703 680 704 CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out 705 706 REAL(wp) :: ze3, zhk 707 REAL(wp), DIMENSION(:,:), POINTER :: zikt 708 709 INTEGER :: ji,jj,jk 710 INTEGER :: ikt,ikb 711 INTEGER, DIMENSION(:,:), POINTER :: mkt, mkb 712 713 CALL wrk_alloc( jpi,jpj, mkt, mkb ) 714 CALL wrk_alloc( jpi,jpj, zikt ) 715 716 ! get first and last level of tbl 717 mkt(:,:) = misfkt(:,:) 718 mkb(:,:) = misfkb(:,:) 719 720 varout(:,:)=0._wp 721 DO jj = 2,jpj 722 DO ji = 2,jpi 723 IF (ssmask(ji,jj) == 1) THEN 724 ikt = mkt(ji,jj) 725 ikb = mkb(ji,jj) 681 ! initialisation 682 pvarout(:,:)=0._wp 683 684 SELECT CASE ( cd_ptin ) 685 CASE ( 'U' ) ! compute U in the top boundary layer at T- point 686 DO jj = 1,jpj 687 DO ji = 1,jpi 688 ikt = miku(ji,jj) ; ikb = miku(ji,jj) 689 ! thickness of boundary layer at least the top level thickness 690 zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3u_n(ji,jj,ikt)) 691 692 ! determine the deepest level influenced by the boundary layer 693 DO jk = ikt+1, mbku(ji,jj) 694 IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 695 END DO 696 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3u_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 697 698 ! level fully include in the ice shelf boundary layer 699 DO jk = ikt, ikb - 1 700 ze3 = fse3u_n(ji,jj,jk) 701 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 702 END DO 703 704 ! level partially include in ice shelf boundary layer 705 zhk = SUM( fse3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 706 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 707 END DO 708 END DO 709 DO jj = 2,jpj 710 DO ji = 2,jpi 711 pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj)) 712 END DO 713 END DO 714 CALL lbc_lnk(pvarout,'T',-1.) 715 716 CASE ( 'V' ) ! compute V in the top boundary layer at T- point 717 DO jj = 1,jpj 718 DO ji = 1,jpi 719 ikt = mikv(ji,jj) ; ikb = mikv(ji,jj) 720 ! thickness of boundary layer at least the top level thickness 721 zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3v_n(ji,jj,ikt)) 722 723 ! determine the deepest level influenced by the boundary layer 724 DO jk = ikt+1, mbkv(ji,jj) 725 IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 726 END DO 727 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3v_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 728 729 ! level fully include in the ice shelf boundary layer 730 DO jk = ikt, ikb - 1 731 ze3 = fse3v_n(ji,jj,jk) 732 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 733 END DO 734 735 ! level partially include in ice shelf boundary layer 736 zhk = SUM( fse3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 737 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 738 END DO 739 END DO 740 DO jj = 2,jpj 741 DO ji = 2,jpi 742 pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1)) 743 END DO 744 END DO 745 CALL lbc_lnk(pvarout,'T',-1.) 746 747 CASE ( 'T' ) ! compute T in the top boundary layer at T- point 748 DO jj = 1,jpj 749 DO ji = 1,jpi 750 ikt = misfkt(ji,jj) 751 ikb = misfkb(ji,jj) 726 752 727 753 ! level fully include in the ice shelf boundary layer 728 754 DO jk = ikt, ikb - 1 729 755 ze3 = fse3t_n(ji,jj,jk) 730 IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 731 IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) & 732 & * r1_hisf_tbl(ji,jj) * ze3 733 IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) & 734 & * r1_hisf_tbl(ji,jj) * ze3 756 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 735 757 END DO 736 758 737 759 ! level partially include in ice shelf boundary layer 738 760 zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 739 IF (cptin == 'T') & 740 & varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk) 741 IF (cptin == 'U') & 742 & varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk) 743 IF (cptin == 'V') & 744 & varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk) 745 END IF 746 END DO 747 END DO 748 749 CALL wrk_dealloc( jpi,jpj, mkt, mkb ) 750 CALL wrk_dealloc( jpi,jpj, zikt ) 751 752 IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.) 753 IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.) 761 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 762 END DO 763 END DO 764 END SELECT 765 766 ! mask mean tbl value 767 pvarout(:,:) = pvarout(:,:) * ssmask(:,:) 768 769 ! deallocation 770 CALL wrk_dealloc( jpi,jpj, zhisf_tbl ) 754 771 ! 755 772 END SUBROUTINE sbc_isf_tbl … … 768 785 !! ** Action : phdivn decreased by the runoff inflow 769 786 !!---------------------------------------------------------------------- 770 REAL(wp), DIMENSION(:,:,:), INTENT( inout) :: phdivn ! horizontal divergence771 ! !787 REAL(wp), DIMENSION(:,:,:), INTENT( inout ) :: phdivn ! horizontal divergence 788 ! 772 789 INTEGER :: ji, jj, jk ! dummy loop indices 773 790 INTEGER :: ikt, ikb 774 INTEGER :: nk_isf 775 REAL(wp) :: zhk, z1_hisf_tbl, zhisf_tbl 776 REAL(wp) :: zfact ! local scalar 791 REAL(wp) :: zhk 792 REAL(wp) :: zfact ! local scalar 777 793 !!---------------------------------------------------------------------- 778 794 ! 779 795 zfact = 0.5_wp 780 796 ! 781 IF ( lk_vvl) THEN ! need to re compute level distribution of isf fresh water797 IF ( lk_vvl ) THEN ! need to re compute level distribution of isf fresh water 782 798 DO jj = 1,jpj 783 799 DO ji = 1,jpi … … 788 804 789 805 ! determine the deepest level influenced by the boundary layer 790 ! test on tmask useless ????? 791 DO jk = ikt, mbkt(ji,jj) 792 IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 806 DO jk = ikt+1, mbkt(ji,jj) 807 IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 793 808 END DO 794 809 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. … … 796 811 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 797 812 798 zhk = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 813 zhk = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 799 814 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer 800 815 END DO 801 816 END DO 802 END IF ! vvl case 803 ! 817 END IF 818 ! 819 !== ice shelf melting distributed over several levels ==! 804 820 DO jj = 1,jpj 805 821 DO ji = 1,jpi … … 809 825 DO jk = ikt, ikb - 1 810 826 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) & 811 & 827 & * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 812 828 END DO 813 829 ! level partially include in ice shelf boundary layer 814 830 phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) & 815 & + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 816 !== ice shelf melting mass distributed over several levels ==! 831 & + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 817 832 END DO 818 833 END DO 819 834 ! 820 835 END SUBROUTINE sbc_isf_div 821 822 823 FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti ) 824 !!---------------------------------------------------------------------- 825 !! *** ROUTINE eos_init *** 826 !! 827 !! ** Purpose : Compute the in-situ temperature [Celcius] 828 !! 829 !! ** Method : 830 !! 831 !! Reference : Bryden,h.,1973,deep-sea res.,20,401-408 832 !!---------------------------------------------------------------------- 833 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptem ! potential temperature [Celcius] 834 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu] 835 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ppress ! pressure [dBar] 836 REAL(wp), DIMENSION(:,:), POINTER :: pti ! in-situ temperature [Celcius] 837 ! REAL(wp) :: fsatg 838 ! REAL(wp) :: pfps, pfpt, pfphp 839 REAL(wp) :: zt, zs, zp, zh, zq, zxk 840 INTEGER :: ji, jj ! dummy loop indices 841 ! 842 CALL wrk_alloc( jpi,jpj, pti ) 843 ! 844 DO jj=1,jpj 845 DO ji=1,jpi 846 zh = ppress(ji,jj) 847 ! Theta1 848 zt = ptem(ji,jj) 849 zs = psal(ji,jj) 850 zp = 0.0 851 zxk= zh * fsatg( zs, zt, zp ) 852 zt = zt + 0.5 * zxk 853 zq = zxk 854 ! Theta2 855 zp = zp + 0.5 * zh 856 zxk= zh*fsatg( zs, zt, zp ) 857 zt = zt + 0.29289322 * ( zxk - zq ) 858 zq = 0.58578644 * zxk + 0.121320344 * zq 859 ! Theta3 860 zxk= zh * fsatg( zs, zt, zp ) 861 zt = zt + 1.707106781 * ( zxk - zq ) 862 zq = 3.414213562 * zxk - 4.121320344 * zq 863 ! Theta4 864 zp = zp + 0.5 * zh 865 zxk= zh * fsatg( zs, zt, zp ) 866 pti(ji,jj) = zt + ( zxk - 2.0 * zq ) / 6.0 867 END DO 868 END DO 869 ! 870 CALL wrk_dealloc( jpi,jpj, pti ) 871 ! 872 END FUNCTION tinsitu 873 874 875 FUNCTION fsatg( pfps, pfpt, pfphp ) 876 !!---------------------------------------------------------------------- 877 !! *** FUNCTION fsatg *** 878 !! 879 !! ** Purpose : Compute the Adiabatic laspse rate [Celcius].[decibar]^-1 880 !! 881 !! ** Reference : Bryden,h.,1973,deep-sea res.,20,401-408 882 !! 883 !! ** units : pressure pfphp decibars 884 !! temperature pfpt deg celsius (ipts-68) 885 !! salinity pfps (ipss-78) 886 !! adiabatic fsatg deg. c/decibar 887 !!---------------------------------------------------------------------- 888 REAL(wp) :: pfps, pfpt, pfphp 889 REAL(wp) :: fsatg 890 ! 891 fsatg = (((-2.1687e-16*pfpt+1.8676e-14)*pfpt-4.6206e-13)*pfphp & 892 & +((2.7759e-12*pfpt-1.1351e-10)*(pfps-35.)+((-5.4481e-14*pfpt & 893 & +8.733e-12)*pfpt-6.7795e-10)*pfpt+1.8741e-8))*pfphp & 894 & +(-4.2393e-8*pfpt+1.8932e-6)*(pfps-35.) & 895 & +((6.6228e-10*pfpt-6.836e-8)*pfpt+8.5258e-6)*pfpt+3.5803e-5 896 ! 897 END FUNCTION fsatg 898 !!====================================================================== 836 !!====================================================================== 899 837 END MODULE sbcisf -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
r6005 r6012 90 90 NAMELIST/namsbc/ nn_fsbc , ln_ana , ln_flx, ln_blk_clio, ln_blk_core, ln_mixcpl, & 91 91 & ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc , ln_rnf , & 92 & ln_ssr , nn_isf , nn_fwb, ln_cdgw , ln_wave , ln_sdw , &92 & ln_ssr , ln_isf , nn_fwb, ln_cdgw , ln_wave , ln_sdw , & 93 93 & nn_lsm , nn_limflx , nn_components, ln_cpl 94 94 INTEGER :: ios … … 143 143 WRITE(numout,*) ' daily mean to diurnal cycle qsr ln_dm2dc = ', ln_dm2dc 144 144 WRITE(numout,*) ' runoff / runoff mouths ln_rnf = ', ln_rnf 145 WRITE(numout,*) ' iceshelf formulation nn_isf = ', nn_isf145 WRITE(numout,*) ' iceshelf melting ln_isf = ', ln_isf 146 146 WRITE(numout,*) ' Sea Surface Restoring on SST and/or SSS ln_ssr = ', ln_ssr 147 147 WRITE(numout,*) ' FreshWater Budget control (=0/1/2) nn_fwb = ', nn_fwb … … 181 181 182 182 ! ! Checks: 183 IF( nn_isf .EQ. 0) THEN ! variable initialisation if no ice shelf184 IF( sbc_isf_alloc() /= 0 ) CALL ctl_stop( ' sbc_init : unable to allocate sbc_isf arrays' )183 IF( .NOT. ln_isf ) THEN ! variable initialisation if no ice shelf 184 IF( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' ) 185 185 fwfisf (:,:) = 0.0_wp ; fwfisf_b (:,:) = 0.0_wp 186 186 risf_tsc(:,:,:) = 0.0_wp ; risf_tsc_b(:,:,:) = 0.0_wp 187 rdivisf = 0.0_wp188 187 END IF 188 189 189 IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa ) fr_i(:,:) = 0.e0 ! no ice in the domain, ice fraction is always zero 190 190 … … 378 378 IF( ln_icebergs ) CALL icb_stp( kt ) ! compute icebergs 379 379 380 IF( nn_isf /= 0 ) CALL sbc_isf( kt )! compute iceshelves380 IF( ln_isf ) CALL sbc_isf( kt ) ! compute iceshelves 381 381 382 382 IF( ln_rnf ) CALL sbc_rnf( kt ) ! add runoffs to fresh water fluxes -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
r5836 r6012 142 142 IF( ln_rnf_tem ) THEN ! use runoffs temperature data 143 143 rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0 144 CALL eos_fzp( sss_m(:,:), ztfrz(:,:) ) 144 145 WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp ) ! if missing data value use SST as runoffs temperature 145 146 rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 146 147 END WHERE 147 148 WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp ) ! where fwf comes from melting of ice shelves or iceberg 148 ztfrz(:,:) = -1.9 !tfreez( sss_m(:,:) ) !PM to be discuss (trouble if sensitivity study) 149 rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * lfusisf * r1_rau0_rcp 149 rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * rlfusisf * r1_rau0_rcp 150 150 END WHERE 151 151 ELSE ! use SST as runoffs temperature -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
r5930 r6012 243 243 ENDIF 244 244 IF( ln_isfcav ) THEN ! ice-shelf cavities 245 IF( ln_traadv_cen .AND. nn_cen_v /= 4 .OR. & ! NO 4th order with ISF246 & ln_traadv_fct .AND. nn_fct_v /= 4 ) CALL ctl_stop( 'tra_adv_init: 4th order COMPACT scheme not allowed with ISF' )245 IF( ln_traadv_cen .AND. nn_cen_v == 4 .OR. & ! NO 4th order with ISF 246 & ln_traadv_fct .AND. nn_fct_v == 4 ) CALL ctl_stop( 'tra_adv_init: 4th order COMPACT scheme not allowed with ISF' ) 247 247 ENDIF 248 248 ! -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
r5836 r6012 103 103 ! 104 104 INTEGER :: ji, jj, jk, jn ! dummy loop indices 105 INTEGER :: ikt 105 106 INTEGER :: ierr ! local integer 106 107 REAL(wp) :: zmsku, zahu_w, zabe1, zcof1, zcoef3 ! local scalars … … 228 229 DO jj = 1, jpjm1 ! bottom correction (partial bottom cell) 229 230 DO ji = 1, fs_jpim1 ! vector opt. 230 !!gm the following anonymous remark is to considered: ! IF useless if zpshde defines pgu everywhere231 231 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 232 232 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) … … 255 255 ELSE ; zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 256 256 ENDIF 257 !!gm I don't understand why we should need this.... since wmask is used instead of tmask258 ! IF ( ln_isfcav ) THEN259 ! DO jj = 1, jpj260 ! DO ji = 1, jpi ! vector opt.261 ! ikt = mikt(ji,jj) ! surface level262 ! zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1)263 ! zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt)264 ! END DO265 ! END DO266 ! END IF267 !!gm268 257 269 258 DO jj = 1 , jpjm1 !== Horizontal fluxes … … 272 261 zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk) 273 262 ! 274 zmsku = 1. / MAX( tmask(ji+1,jj,jk ) + tmask(ji,jj,jk+1) &275 & + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk ), 1. )276 ! 277 zmskv = 1. / MAX( tmask(ji,jj+1,jk ) + tmask(ji,jj,jk+1) &278 & + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk ), 1. )263 zmsku = 1. / MAX( wmask(ji+1,jj,jk ) + wmask(ji,jj,jk+1) & 264 & + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk ), 1. ) 265 ! 266 zmskv = 1. / MAX( wmask(ji,jj+1,jk ) + wmask(ji,jj,jk+1) & 267 & + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk ), 1. ) 279 268 ! 280 269 zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku … … 317 306 DO ji = fs_2, fs_jpim1 ! vector opt. 318 307 ! 319 zmsku = tmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) &308 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 320 309 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 321 zmskv = tmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) &310 zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 322 311 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 323 312 ! … … 343 332 DO jj = 1, jpjm1 344 333 DO ji = fs_2, fs_jpim1 345 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk) &334 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * wmask(ji,jj,jk) & 346 335 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & 347 336 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) … … 358 347 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) & 359 348 & + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj) & 360 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) / fse3w(ji,jj,jk)349 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * wmask(ji,jj,jk) / fse3w(ji,jj,jk) 361 350 END DO 362 351 END DO … … 366 355 DO jj = 1, jpjm1 367 356 DO ji = fs_2, fs_jpim1 368 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk) &357 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * wmask(ji,jj,jk) & 369 358 & * ( ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) ) & 370 359 & + akz (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) ) ) -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
r5930 r6012 28 28 USE sbc_oce ! surface boundary condition: ocean 29 29 USE sbcrnf ! river runoffs 30 USE sbcisf ! ice shelf melting /freezing30 USE sbcisf ! ice shelf melting 31 31 USE zdf_oce ! ocean vertical mixing 32 32 USE domvvl ! variable volume … … 262 262 ll_traqsr = ln_traqsr ! active tracers case and solar penetration 263 263 ll_rnf = ln_rnf ! active tracers case and river runoffs 264 IF (nn_isf .GE. 1) THEN 265 ll_isf = .TRUE. ! active tracers case and ice shelf melting/freezing 266 ELSE 267 ll_isf = .FALSE. 268 END IF 264 ll_isf = ln_isf ! active tracers case and ice shelf melting 269 265 ELSE 270 266 ll_traqsr = .FALSE. ! active tracers case and NO solar penetration 271 267 ll_rnf = .FALSE. ! passive tracers or NO river runoffs 272 ll_isf = .FALSE. ! passive tracers or NO ice shelf melting /freezing268 ll_isf = .FALSE. ! passive tracers or NO ice shelf melting 273 269 ENDIF 274 270 ! -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
r6006 r6012 119 119 INTEGER :: ikt, ikb 120 120 REAL(wp) :: zfact, z1_e3t, zdep 121 REAL(wp) :: z alpha, zhk121 REAL(wp) :: zt_frz, zpress 122 122 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds 123 123 !!---------------------------------------------------------------------- … … 218 218 !---------------------------------------- 219 219 ! 220 IF( nn_isf > 0) THEN220 IF( ln_isf ) THEN 221 221 zfact = 0.5_wp 222 222 DO jj = 2, jpj … … 227 227 228 228 ! level fully include in the ice shelf boundary layer 229 ! if isfdiv, we have to remove heat flux due to inflow at 0oC (as in rnf when you add rnf at sst)230 229 ! sign - because fwf sign of evapo (rnf sign of precip) 231 230 DO jk = ikt, ikb - 1 232 231 ! compute trend 233 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & 234 & + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) 235 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) & 236 & + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) 232 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & 233 & + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) ) & 234 & * r1_hisf_tbl(ji,jj) 237 235 END DO 238 236 239 237 ! level partially include in ice shelf boundary layer 240 238 ! compute trend 241 tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem) &242 & + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)243 tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal) &244 & + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 239 tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem) & 240 & + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) ) & 241 & * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 242 245 243 END DO 246 244 END DO -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
r5836 r6012 191 191 ! 192 192 END SUBROUTINE zps_hde 193 194 195 SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu , pgtv , pgtui, pgtvi, & 196 & prd, pgru , pgrv , pmru , pmrv , pgzu , pgzv , pge3ru , pge3rv , & 197 & pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi ) 198 !!---------------------------------------------------------------------- 199 !! *** ROUTINE zps_hde *** 193 ! 194 SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, pgtui, pgtvi, & 195 & prd, pgru, pgrv, pgrui, pgrvi ) 196 !!---------------------------------------------------------------------- 197 !! *** ROUTINE zps_hde_isf *** 200 198 !! 201 199 !! ** Purpose : Compute the horizontal derivative of T, S and rho 202 200 !! at u- and v-points with a linear interpolation for z-coordinate 203 !! with partial steps .201 !! with partial steps for top (ice shelf) and bottom. 204 202 !! 205 203 !! ** Method : In z-coord with partial steps, scale factors on last 206 204 !! levels are different for each grid point, so that T, S and rd 207 205 !! points are not at the same depth as in z-coord. To have horizontal 208 !! gradients again, we interpolate T and S at the good depth : 206 !! gradients again, we interpolate T and S at the good depth : 207 !! For the bottom case: 209 208 !! Linear interpolation of T, S 210 209 !! Computation of di(tb) and dj(tb) by vertical interpolation: … … 235 234 !! di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~ 236 235 !! 236 !! For the top case (ice shelf): As for the bottom case but upside down 237 !! 237 238 !! ** Action : compute for top and bottom interfaces 238 239 !! - pgtu, pgtv, pgtui, pgtvi: horizontal gradient of tracer at u- & v-points 239 240 !! - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points 240 !! - pmru, pmrv, pmrui, pmrvi: horizontal sum of rho at u- & v- point (used in dynhpg with vvl) 241 !! - pgzu, pgzv, pgzui, pgzvi: horizontal gradient of z at u- and v- point (used in dynhpg with vvl) 242 !! - pge3ru, pge3rv, pge3rui, pge3rvi: horizontal gradient of rho weighted by local e3w at u- & v-points 243 !!---------------------------------------------------------------------- 244 INTEGER , INTENT(in ) :: kt ! ocean time-step index 245 INTEGER , INTENT(in ) :: kjpt ! number of tracers 246 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pta ! 4D tracers fields 247 ! !! u-point ! v-point ! 248 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtu , pgtv ! bottom GRADh( ptra ) 249 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtui , pgtvi ! top GRADh( ptra ) 250 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 251 ! !! u-point ! v-point ! 252 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru , pgrv ! bottom GRADh( prd ) 253 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pmru , pmrv ! bottom SUM ( prd ) 254 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgzu , pgzv ! bottom GRADh( z ) 255 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pge3ru , pge3rv ! bottom GRADh( prd ) weighted by e3w 256 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgrui , pgrvi ! top GRADh( prd ) 257 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pmrui , pmrvi ! top SUM ( prd ) 258 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgzui , pgzvi ! top GRADh( z ) 259 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pge3rui , pge3rvi ! top GRADh( prd ) weighted by e3w 241 !!---------------------------------------------------------------------- 242 INTEGER , INTENT(in ) :: kt ! ocean time-step index 243 INTEGER , INTENT(in ) :: kjpt ! number of tracers 244 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pta ! 4D tracers fields 245 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 246 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF) 247 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 248 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 249 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top) 260 250 ! 261 251 INTEGER :: ji, jj, jn ! Dummy loop indices 262 252 INTEGER :: iku, ikv, ikum1, ikvm1,ikup1, ikvp1 ! partial step level (ocean bottom level) at u- and v-points 263 REAL(wp) :: ze3wu, ze3wv, zmaxu, zmaxv , zdzwu, zdzwv, zdzwuip1, zdzwvjp1! temporary scalars253 REAL(wp) :: ze3wu, ze3wv, zmaxu, zmaxv ! temporary scalars 264 254 REAL(wp), DIMENSION(jpi,jpj) :: zri, zrj, zhi, zhj ! NB: 3rd dim=1 to use eos 265 255 REAL(wp), DIMENSION(jpi,jpj,kjpt) :: zti, ztj ! … … 277 267 DO jj = 1, jpjm1 278 268 DO ji = 1, jpim1 279 iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 280 ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 ) ! if level first is a p-step, ik.m1=1 269 270 iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 271 ikv = mbkv(ji,jj); ikvm1 = MAX( ikv - 1 , 1 ) ! if level first is a p-step, ik.m1=1 272 ze3wu = fsdept_n(ji+1,jj,iku) - fsdept_n(ji,jj,iku) 273 ze3wv = fsdept_n(ji,jj+1,ikv) - fsdept_n(ji,jj,ikv) 274 ! 275 ! i- direction 276 IF( ze3wu >= 0._wp ) THEN ! case 1 277 zmaxu = ze3wu / fse3w(ji+1,jj,iku) 278 ! interpolated values of tracers 279 zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 280 ! gradient of tracers 281 pgtu(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 282 ELSE ! case 2 283 zmaxu = -ze3wu / fse3w(ji,jj,iku) 284 ! interpolated values of tracers 285 zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 286 ! gradient of tracers 287 pgtu(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 288 ENDIF 289 ! 290 ! j- direction 291 IF( ze3wv >= 0._wp ) THEN ! case 1 292 zmaxv = ze3wv / fse3w(ji,jj+1,ikv) 293 ! interpolated values of tracers 294 ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 295 ! gradient of tracers 296 pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 297 ELSE ! case 2 298 zmaxv = -ze3wv / fse3w(ji,jj,ikv) 299 ! interpolated values of tracers 300 ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 301 ! gradient of tracers 302 pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 303 ENDIF 304 305 END DO 306 END DO 307 CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. ) ! Lateral boundary cond. 308 ! 309 END DO 310 311 ! horizontal derivative of density anomalies (rd) 312 IF( PRESENT( prd ) ) THEN ! depth of the partial step level 313 pgru(:,:)=0.0_wp ; pgrv(:,:)=0.0_wp ; 314 ! 315 DO jj = 1, jpjm1 316 DO ji = 1, jpim1 317 318 iku = mbku(ji,jj) 319 ikv = mbkv(ji,jj) 320 ze3wu = fsdept_n(ji+1,jj,iku) - fsdept_n(ji,jj,iku) 321 ze3wv = fsdept_n(ji,jj+1,ikv) - fsdept_n(ji,jj,ikv) 322 ! 323 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = fsdept(ji ,jj,iku) ! i-direction: case 1 324 ELSE ; zhi(ji,jj) = fsdept(ji+1,jj,iku) ! - - case 2 325 ENDIF 326 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = fsdept(ji,jj ,ikv) ! j-direction: case 1 327 ELSE ; zhj(ji,jj) = fsdept(ji,jj+1,ikv) ! - - case 2 328 ENDIF 329 330 END DO 331 END DO 332 333 ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 334 ! step and store it in zri, zrj for each case 335 CALL eos( zti, zhi, zri ) 336 CALL eos( ztj, zhj, zrj ) 337 338 DO jj = 1, jpjm1 ! Gradient of density at the last level 339 DO ji = 1, jpim1 340 341 iku = mbku(ji,jj) 342 ikv = mbkv(ji,jj) 343 ze3wu = fsdept_n(ji+1,jj,iku) - fsdept_n(ji,jj,iku) 344 ze3wv = fsdept_n(ji,jj+1,ikv) - fsdept_n(ji,jj,ikv) 345 346 IF( ze3wu >= 0._wp ) THEN ; pgru(ji,jj) = ssumask(ji,jj) * ( zri(ji ,jj ) - prd(ji,jj,iku) ) ! i: 1 347 ELSE ; pgru(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj,iku) - zri(ji,jj ) ) ! i: 2 348 ENDIF 349 IF( ze3wv >= 0._wp ) THEN ; pgrv(ji,jj) = ssvmask(ji,jj) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 350 ELSE ; pgrv(ji,jj) = ssvmask(ji,jj) * ( prd(ji,jj+1,ikv) - zrj(ji,jj ) ) ! j: 2 351 ENDIF 352 353 END DO 354 END DO 355 356 CALL lbc_lnk( pgru , 'U', -1. ) ; CALL lbc_lnk( pgrv , 'V', -1. ) ! Lateral boundary conditions 357 ! 358 END IF 359 ! 360 ! !== (ISH) compute grui and gruvi ==! 361 ! 362 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! ! 363 DO jj = 1, jpjm1 364 DO ji = 1, jpim1 365 iku = miku(ji,jj); ikup1 = miku(ji,jj) + 1 366 ikv = mikv(ji,jj); ikvp1 = mikv(ji,jj) + 1 367 ! 281 368 ! (ISF) case partial step top and bottom in adjacent cell in vertical 282 369 ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 283 370 ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 284 371 ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 285 ze3wu = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))286 ze3wv = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))287 ! 372 ze3wu = fsdept_n(ji,jj,iku) - fsdept_n(ji+1,jj,iku) 373 ze3wv = fsdept_n(ji,jj,ikv) - fsdept_n(ji,jj+1,ikv) 374 288 375 ! i- direction 289 376 IF( ze3wu >= 0._wp ) THEN ! case 1 290 zmaxu = ze3wu / fse3w(ji+1,jj,iku) 291 ! interpolated values of tracers 292 zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 377 zmaxu = ze3wu / fse3w(ji+1,jj,ikup1) 378 ! interpolated values of tracers 379 zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikup1,jn) - pta(ji+1,jj,iku,jn) ) 380 ! gradient of tracers 381 pgtui(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 382 ELSE ! case 2 383 zmaxu = - ze3wu / fse3w(ji,jj,ikup1) 384 ! interpolated values of tracers 385 zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikup1,jn) - pta(ji,jj,iku,jn) ) 293 386 ! gradient of tracers 294 pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 295 ELSE ! case 2 296 zmaxu = -ze3wu / fse3w(ji,jj,iku) 297 ! interpolated values of tracers 298 zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 299 ! gradient of tracers 300 pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 387 pgtui(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 301 388 ENDIF 302 389 ! 303 390 ! j- direction 304 391 IF( ze3wv >= 0._wp ) THEN ! case 1 305 zmaxv = ze3wv / fse3w(ji,jj+1,ikv) 306 ! interpolated values of tracers 307 ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 308 ! gradient of tracers 309 pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 310 ELSE ! case 2 311 zmaxv = -ze3wv / fse3w(ji,jj,ikv) 312 ! interpolated values of tracers 313 ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 314 ! gradient of tracers 315 pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 316 ENDIF 317 END DO 318 END DO 319 CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. ) ! Lateral boundary cond. 392 zmaxv = ze3wv / fse3w(ji,jj+1,ikvp1) 393 ! interpolated values of tracers 394 ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvp1,jn) - pta(ji,jj+1,ikv,jn) ) 395 ! gradient of tracers 396 pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 397 ELSE ! case 2 398 zmaxv = - ze3wv / fse3w(ji,jj,ikvp1) 399 ! interpolated values of tracers 400 ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvp1,jn) - pta(ji,jj,ikv,jn) ) 401 ! gradient of tracers 402 pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 403 ENDIF 404 405 END DO 406 END DO 407 CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. ); CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. ) ! Lateral boundary cond. 320 408 ! 321 409 END DO … … 323 411 IF( PRESENT( prd ) ) THEN !== horizontal derivative of density anomalies (rd) ==! (optional part) 324 412 ! 325 pgru (:,:)=0._wp ; pgrv (:,:) = 0._wp 326 pgzu (:,:)=0._wp ; pgzv (:,:) = 0._wp 327 pmru (:,:)=0._wp ; pmru (:,:) = 0._wp 328 pge3ru(:,:)=0._wp ; pge3rv(:,:) = 0._wp 329 ! 330 DO jj = 1, jpjm1 ! depth of the partial step level 331 DO ji = 1, jpim1 332 iku = mbku(ji,jj) 333 ikv = mbkv(ji,jj) 334 ze3wu = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 335 ze3wv = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 336 ! 337 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu ! i-direction: case 1 338 ELSE ; zhi(ji,jj) = fsdept(ji ,jj,iku) + ze3wu ! - - case 2 339 ENDIF 340 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = fsdept(ji,jj+1,ikv) - ze3wv ! j-direction: case 1 341 ELSE ; zhj(ji,jj) = fsdept(ji,jj ,ikv) + ze3wv ! - - case 2 342 ENDIF 413 pgrui(:,:) =0.0_wp; pgrvi(:,:) =0.0_wp; 414 DO jj = 1, jpjm1 415 DO ji = 1, jpim1 416 417 iku = miku(ji,jj) 418 ikv = mikv(ji,jj) 419 ze3wu = fsdept_n(ji,jj,iku) - fsdept_n(ji+1,jj,iku) 420 ze3wv = fsdept_n(ji,jj,ikv) - fsdept_n(ji,jj+1,ikv) 421 ! 422 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = fsdept(ji ,jj,iku) ! i-direction: case 1 423 ELSE ; zhi(ji,jj) = fsdept(ji+1,jj,iku) ! - - case 2 424 ENDIF 425 426 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = fsdept(ji,jj ,ikv) ! j-direction: case 1 427 ELSE ; zhj(ji,jj) = fsdept(ji,jj+1,ikv) ! - - case 2 428 ENDIF 429 343 430 END DO 344 431 END DO … … 346 433 CALL eos( zti, zhi, zri ) ! interpolated density from zti, ztj 347 434 CALL eos( ztj, zhj, zrj ) ! at the partial step depth output in zri, zrj 348 435 ! 349 436 DO jj = 1, jpjm1 ! Gradient of density at the last level 350 437 DO ji = 1, jpim1 351 iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 352 ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 ) ! last and before last ocean level at u- & v-points 353 ze3wu = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 354 ze3wv = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 355 IF( ze3wu >= 0._wp ) THEN 356 pgzu(ji,jj) = (fsde3w(ji+1,jj,iku) - ze3wu) - fsde3w(ji,jj,iku) 357 pgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji ,jj) - prd(ji,jj,iku) ) ! i: 1 358 pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji ,jj) + prd(ji,jj,iku) ) ! i: 1 359 pge3ru(ji,jj) = umask(ji,jj,iku) & 360 * ( (fse3w(ji+1,jj,iku) - ze3wu )* ( zri(ji ,jj ) + prd(ji+1,jj,ikum1) + 2._wp) & 361 - fse3w(ji ,jj,iku) * ( prd(ji ,jj,iku) + prd(ji ,jj,ikum1) + 2._wp) ) ! j: 2 362 ELSE 363 pgzu(ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) + ze3wu) 364 pgru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) ) ! i: 2 365 pmru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) ) ! i: 2 366 pge3ru(ji,jj) = umask(ji,jj,iku) & 367 * ( fse3w(ji+1,jj,iku) * ( prd(ji+1,jj,iku) + prd(ji+1,jj,ikum1) + 2._wp) & 368 -(fse3w(ji ,jj,iku) + ze3wu) * ( zri(ji ,jj ) + prd(ji ,jj,ikum1) + 2._wp) ) ! j: 2 369 ENDIF 370 IF( ze3wv >= 0._wp ) THEN 371 pgzv(ji,jj) = (fsde3w(ji,jj+1,ikv) - ze3wv) - fsde3w(ji,jj,ikv) 372 pgrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 373 pmrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) + prd(ji,jj,ikv) ) ! j: 1 374 pge3rv(ji,jj) = vmask(ji,jj,ikv) & 375 * ( (fse3w(ji,jj+1,ikv) - ze3wv )* ( zrj(ji,jj ) + prd(ji,jj+1,ikvm1) + 2._wp) & 376 - fse3w(ji,jj ,ikv) * ( prd(ji,jj ,ikv) + prd(ji,jj ,ikvm1) + 2._wp) ) ! j: 2 377 ELSE 378 pgzv(ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) + ze3wv) 379 pgrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) ) ! j: 2 380 pmrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) ) ! j: 2 381 pge3rv(ji,jj) = vmask(ji,jj,ikv) & 382 * ( fse3w(ji,jj+1,ikv) * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikvm1) + 2._wp) & 383 -(fse3w(ji,jj ,ikv) + ze3wv) * ( zrj(ji,jj ) + prd(ji,jj ,ikvm1) + 2._wp) ) ! j: 2 384 ENDIF 385 END DO 386 END DO 387 CALL lbc_lnk( pgru , 'U', -1. ) ; CALL lbc_lnk( pgrv , 'V', -1. ) ! Lateral boundary conditions 388 CALL lbc_lnk( pmru , 'U', 1. ) ; CALL lbc_lnk( pmrv , 'V', 1. ) ! Lateral boundary conditions 389 CALL lbc_lnk( pgzu , 'U', -1. ) ; CALL lbc_lnk( pgzv , 'V', -1. ) ! Lateral boundary conditions 390 CALL lbc_lnk( pge3ru , 'U', -1. ) ; CALL lbc_lnk( pge3rv , 'V', -1. ) ! Lateral boundary conditions 391 ! 392 END IF 393 ! 394 ! !== (ISH) compute grui and gruvi ==! 395 ! 396 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! ! 397 DO jj = 1, jpjm1 398 DO ji = 1, jpim1 399 iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1 400 ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1 401 ! 402 ! (ISF) case partial step top and bottom in adjacent cell in vertical 403 ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 404 ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 405 ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 406 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 407 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 408 ! i- direction 409 IF( ze3wu >= 0._wp ) THEN ! case 1 410 zmaxu = ze3wu / fse3w(ji+1,jj,iku+1) 411 ! interpolated values of tracers 412 zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) ) 413 ! gradient of tracers 414 pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 415 ELSE ! case 2 416 zmaxu = - ze3wu / fse3w(ji,jj,iku+1) 417 ! interpolated values of tracers 418 zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) ) 419 ! gradient of tracers 420 pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 421 ENDIF 422 ! 423 ! j- direction 424 IF( ze3wv >= 0._wp ) THEN ! case 1 425 zmaxv = ze3wv / fse3w(ji,jj+1,ikv+1) 426 ! interpolated values of tracers 427 ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) ) 428 ! gradient of tracers 429 pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 430 ELSE ! case 2 431 zmaxv = - ze3wv / fse3w(ji,jj,ikv+1) 432 ! interpolated values of tracers 433 ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) ) 434 ! gradient of tracers 435 pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 436 ENDIF 437 END DO!! 438 END DO!! 439 CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. ) ! Lateral boundary cond. 440 ! 441 END DO 442 443 IF( PRESENT( prd ) ) THEN !== horizontal derivative of density anomalies (rd) ==! (optional part) 444 ! 445 pgrui(:,:) =0.0_wp ; pgrvi(:,:) =0.0_wp ; 446 pgzui(:,:) =0.0_wp ; pgzvi(:,:) =0.0_wp ; 447 pmrui(:,:) =0.0_wp ; pmrui(:,:) =0.0_wp ; 448 pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ; 449 ! 450 DO jj = 1, jpjm1 ! depth of the partial step level 451 DO ji = 1, jpim1 452 iku = miku(ji,jj) 453 ikv = mikv(ji,jj) 454 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 455 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 456 ! 457 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu ! i-direction: case 1 458 ELSE ; zhi(ji,jj) = fsdept(ji ,jj,iku) - ze3wu ! - - case 2 459 ENDIF 460 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv ! j-direction: case 1 461 ELSE ; zhj(ji,jj) = fsdept(ji,jj ,ikv) - ze3wv ! - - case 2 462 ENDIF 463 END DO 464 END DO 465 ! 466 CALL eos( zti, zhi, zri ) ! interpolated density from zti, ztj 467 CALL eos( ztj, zhj, zrj ) ! at the partial step depth output in zri, zrj 468 ! 469 DO jj = 1, jpjm1 ! Gradient of density at the last level 470 DO ji = 1, jpim1 471 iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1 472 ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1 473 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 474 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 475 IF( ze3wu >= 0._wp ) THEN 476 pgzui (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku) 477 pgrui (ji,jj) = umask(ji,jj,iku) * ( zri(ji,jj) - prd(ji,jj,iku) ) ! i: 1 478 pmrui (ji,jj) = umask(ji,jj,iku) * ( zri(ji,jj) + prd(ji,jj,iku) ) ! i: 1 479 pge3rui(ji,jj) = umask(ji,jj,iku+1) & 480 & * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj ) + prd(ji+1,jj,iku+1) + 2._wp) & 481 & - fse3w(ji ,jj,iku+1) * (prd(ji,jj,iku) + prd(ji ,jj,iku+1) + 2._wp) ) ! i: 1 482 ELSE 483 pgzui (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu) 484 pgrui (ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) ) ! i: 2 485 pmrui (ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) ) ! i: 2 486 pge3rui(ji,jj) = umask(ji,jj,iku+1) & 487 & * ( fse3w(ji+1,jj,iku+1) * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp) & 488 & -(fse3w(ji ,jj,iku+1) + ze3wu) * (zri(ji,jj ) + prd(ji ,jj,iku+1) + 2._wp) ) ! i: 2 489 ENDIF 490 IF( ze3wv >= 0._wp ) THEN 491 pgzvi (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv) 492 pgrvi (ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 493 pmrvi (ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) + prd(ji,jj,ikv) ) ! j: 1 494 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1) & 495 & * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj ) + prd(ji,jj+1,ikv+1) + 2._wp) & 496 & - fse3w(ji,jj ,ikv+1) * ( prd(ji,jj,ikv) + prd(ji,jj ,ikv+1) + 2._wp) ) ! j: 1 497 ! + 2 due to the formulation in density and not in anomalie in hpg sco 498 ELSE 499 pgzvi (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv) 500 pgrvi (ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) ) ! j: 2 501 pmrvi (ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) ) ! j: 2 502 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1) & 503 & * ( fse3w(ji,jj+1,ikv+1) * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) & 504 & -(fse3w(ji,jj ,ikv+1) + ze3wv) * ( zrj(ji,jj ) + prd(ji,jj ,ikv+1) + 2._wp) ) ! j: 2 505 ENDIF 506 END DO 507 END DO 508 CALL lbc_lnk( pgrui , 'U', -1. ) ; CALL lbc_lnk( pgrvi , 'V', -1. ) ! Lateral boundary conditions 509 CALL lbc_lnk( pmrui , 'U', 1. ) ; CALL lbc_lnk( pmrvi , 'V', 1. ) ! Lateral boundary conditions 510 CALL lbc_lnk( pgzui , 'U', -1. ) ; CALL lbc_lnk( pgzvi , 'V', -1. ) ! Lateral boundary conditions 511 CALL lbc_lnk( pge3rui , 'U', -1. ) ; CALL lbc_lnk( pge3rvi , 'V', -1. ) ! Lateral boundary conditions 438 439 iku = miku(ji,jj) 440 ikv = mikv(ji,jj) 441 ze3wu = fsdept_n(ji,jj,iku) - fsdept_n(ji+1,jj,iku) 442 ze3wv = fsdept_n(ji,jj,ikv) - fsdept_n(ji,jj+1,ikv) 443 444 IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji ,jj ) - prd(ji,jj,iku) ) ! i: 1 445 ELSE ; pgrui(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj ,iku) - zri(ji,jj ) ) ! i: 2 446 ENDIF 447 448 IF( ze3wv >= 0._wp ) THEN ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( zrj(ji ,jj ) - prd(ji,jj,ikv) ) ! j: 1 449 ELSE ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( prd(ji ,jj+1,ikv) - zrj(ji,jj ) ) ! j: 2 450 ENDIF 451 452 END DO 453 END DO 454 CALL lbc_lnk( pgrui , 'U', -1. ); CALL lbc_lnk( pgrvi , 'V', -1. ) ! Lateral boundary conditions 512 455 ! 513 456 END IF -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
r4990 r6012 31 31 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmld !: mixing layer depth (turbocline) [m] 32 32 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlp !: mixed layer depth (rho=rho0+zdcrit) [m] 33 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlpt !: mixed layer depth at t-points [m]33 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlpt !: depth of the last T-point inside the mixed layer 34 34 35 35 REAL(wp), PUBLIC :: rho_c = 0.01_wp !: density criterion for mixed layer depth … … 111 111 END DO 112 112 ! 113 ! w-level of the turbocline 113 ! w-level of the turbocline and mixing layer (iom_use) 114 114 imld(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point 115 115 DO jk = jpkm1, nlb10, -1 ! from the bottom to nlb10 116 116 DO jj = 1, jpj 117 117 DO ji = 1, jpi 118 imkt = mikt(ji,jj) 119 IF( avt (ji,jj,jk) < avt_c ) imld(ji,jj) = MAX( imkt, jk ) ! Turbocline 118 IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) ) imld(ji,jj) = jk ! Turbocline 120 119 END DO 121 120 END DO … … 127 126 iikn = nmln(ji,jj) 128 127 imkt = mikt(ji,jj) 129 hmld (ji,jj) = ( fsdepw(ji,jj,iiki ) - fsdepw(ji,jj,imkt ) )* ssmask(ji,jj) ! Turbocline depth130 hmlp (ji,jj) = ( fsdepw(ji,jj,iikn ) - fsdepw(ji,jj,MAX( imkt,nla10 ) )) * ssmask(ji,jj) ! Mixed layer depth131 hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt ) )* ssmask(ji,jj) ! depth of the last T-point inside the mixed layer128 hmld (ji,jj) = fsdepw(ji,jj,iiki ) * ssmask(ji,jj) ! Turbocline depth 129 hmlp (ji,jj) = fsdepw(ji,jj,iikn ) * ssmask(ji,jj) ! Mixed layer depth 130 hmlpt(ji,jj) = fsdept(ji,jj,iikn-1) * ssmask(ji,jj) ! depth of the last T-point inside the mixed layer 132 131 END DO 133 132 END DO 134 133 IF( .NOT.lk_offline ) THEN ! no need to output in offline mode 135 CALL iom_put( "mldr10_1", hmlp ) ! mixed layer depth 136 CALL iom_put( "mldkz5" , hmld ) ! turbocline depth 134 IF ( iom_use("mldr10_1") ) THEN 135 IF( .NOT. ln_isfcav ) CALL iom_put( "mldr10_1", hmlp ) ! mixed layer depth 136 IF( ln_isfcav ) CALL iom_put( "mldr10_1", hmlp - risfdep) ! mixed layer thickness 137 END IF 138 IF ( iom_use("mldkz5") ) THEN 139 IF( .NOT. ln_isfcav ) CALL iom_put( "mldkz5" , hmld ) ! turbocline depth 140 IF( ln_isfcav ) CALL iom_put( "mldkz5" , hmld - risfdep ) ! turbocline thickness 141 END IF 137 142 ENDIF 138 143 -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/oce.F90
r5930 r6012 60 60 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gtsu, gtsv !: horizontal gradient of T, S bottom u-point 61 61 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gru , grv !: horizontal gradient of rd at bottom u-point 62 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: aru , arv63 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gzu , gzv64 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ge3ru, ge3rv !: horizontal gradient of T, S and rd at top v-point65 62 66 63 !! (ISF) interpolated gradient (only used for ice shelf case) … … 68 65 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gtui, gtvi !: horizontal gradient of T, S and rd at top u-point 69 66 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: grui, grvi !: horizontal gradient of T, S and rd at top v-point 70 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: arui, arvi !: horizontal average of rd at top v-point71 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gzui, gzvi !: horizontal gradient of z at top v-point72 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ge3rui, ge3rvi !: horizontal gradient of T, S and rd at top v-point73 67 !! (ISF) ice load 74 68 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: riceload … … 115 109 & spgu (jpi,jpj) , spgv(jpi,jpj) , & 116 110 & gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts), & 117 & aru(jpi,jpj) , arv(jpi,jpj) , &118 & gzu(jpi,jpj) , gzv(jpi,jpj) , &119 111 & gru(jpi,jpj) , grv(jpi,jpj) , & 120 & ge3ru(jpi,jpj) , ge3rv(jpi,jpj) , &121 112 & gtui(jpi,jpj,jpts), gtvi(jpi,jpj,jpts), & 122 & arui(jpi,jpj) , arvi(jpi,jpj) , &123 & gzui(jpi,jpj) , gzvi(jpi,jpj) , &124 & ge3rui(jpi,jpj) , ge3rvi(jpi,jpj) , &125 113 & grui(jpi,jpj) , grvi(jpi,jpj) , & 126 114 & riceload(jpi,jpj), STAT=ierr(2) ) -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/step.F90
r6010 r6012 166 166 167 167 IF( ln_zps .AND. ln_isfcav) & 168 & CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 169 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 170 & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the first ocean level 171 168 & CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 169 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level 172 170 IF( ln_traldf_triad ) THEN 173 171 CALL ldf_slp_triad( kstp ) ! before slope for triad operator … … 186 184 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 187 185 CALL wzv ( kstp ) ! now cross-level velocity 188 189 186 !!gm : why also here ???? 190 187 IF(ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations … … 196 193 !! but ensures reproductible results 197 194 !! with previous versions using split-explicit free surface 198 IF( ln_zps .AND. .NOT. ln_isfcav) & ! Partial steps: bottom before horizontal gradient 199 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! of t, s, rd at the last ocean level 200 & rhd, gru , grv ) 201 IF( ln_zps .AND. ln_isfcav) & ! Partial steps: top & bottom before horizontal gradient 202 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi, & 203 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 204 & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) 195 IF( ln_zps .AND. .NOT. ln_isfcav ) & 196 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient 197 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 198 IF( ln_zps .AND. ln_isfcav ) & 199 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 200 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level 205 201 !!jc: fs simplification 206 202 -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/TOP_SRC/trcini.F90
r5836 r6012 246 246 END SUBROUTINE trc_ini_state 247 247 248 249 248 SUBROUTINE top_alloc 250 249 !!---------------------------------------------------------------------- -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/SETTE/input_ORCA2_LIM.cfg
r5398 r6012 1 ORCA2_LIM_nemo_v3.6 .tar ORCA2_LIM_nemo_v3.61 ORCA2_LIM_nemo_v3.6st.tar ORCA2_LIM_nemo_v3.6
Note: See TracChangeset
for help on using the changeset viewer.