Changeset 5098 for branches/2015/dev_r5094_UKMO_ISFCLEAN
- Timestamp:
- 2015-02-20T20:11:47+01:00 (9 years ago)
- Location:
- branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO
- Files:
-
- 33 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90
r4990 r5098 537 537 CALL eos_rab( pts, rab_n ) ! now local thermal/haline expension ratio at T-points 538 538 CALL bn2 ( pts, rab_n, rn2 ) ! now Brunt-Vaisala 539 IF( ln_zps ) & ! Partial steps: before Horizontal DErivative 540 & CALL zps_hde( kt, jpts, pts, gtsu, gtsv, & ! Partial steps: before horizontal gradient 541 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & ! 542 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 539 ! Partial steps: before Horizontal DErivative 543 540 ! only gtsu, gtsv, rhd, gru , grv are used 541 IF( ln_zps ) CALL zps_hde ( kt, jpts, pts, gtsu, gtsv, & ! Partial steps: before horizontal gradient 542 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 543 IF( ln_zps .AND. ln_isfcav) & 544 & CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, & ! Partial steps for top cell (ISF) 545 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 546 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the first ocean level 547 544 548 545 549 -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r4998 r5098 746 746 747 747 748 IF( ln_zps .AND. .NOT. lk_c1d ) & 749 & CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 750 & rhd, gru , grv, aru, arv, gzu, gzv, ge3ru, ge3rv, & ! 751 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 748 IF( ln_zps .AND. .NOT. lk_c1d ) & 749 & CALL zps_hde ( kt, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 750 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 751 IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) & 752 & CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps for top cell (ISF) 753 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 754 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 752 755 753 756 #if defined key_zdfkpp -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
r4990 r5098 145 145 z3d(:,:,:) = tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) 146 146 CALL iom_put( "toce" , z3d ) ! heat content 147 DO jj = 1, jpj 148 DO ji = 1, jpi 149 z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * fse3t_n(ji,jj,mikt(ji,jj)) 150 END DO 151 END DO 152 CALL iom_put( "sst" , z2d(:,:) ) ! sea surface heat content 153 DO jj = 1, jpj 154 DO ji = 1, jpi 155 z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem)**2 * fse3t_n(ji,jj,mikt(ji,jj)) 156 END DO 157 END DO 158 CALL iom_put( "sst2" , z2d(:,:) ) ! sea surface content of squared temperature 147 CALL iom_put( "sst" , z3d(:,:,1) ) ! sea surface heat content 148 z3d(:,:,1) = tsn(:,:,1,jp_tem) * z3d(:,:,1) 149 CALL iom_put( "sst2" , z3d(:,:,1) ) ! sea surface content of squared temperature 159 150 z3d(:,:,:) = tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) 160 151 CALL iom_put( "soce" , z3d ) ! salinity content 161 DO jj = 1, jpj 162 DO ji = 1, jpi 163 z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * fse3t_n(ji,jj,mikt(ji,jj)) 164 END DO 165 END DO 166 CALL iom_put( "sss" , z2d(:,:) ) ! sea surface salinity content 167 DO jj = 1, jpj 168 DO ji = 1, jpi 169 z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal)**2 * fse3t_n(ji,jj,mikt(ji,jj)) 170 END DO 171 END DO 172 CALL iom_put( "sss2" , z2d(:,:) ) ! sea surface content of squared salinity 152 CALL iom_put( "sss" , z3d(:,:,1) ) ! sea surface salinity content 153 z3d(:,:,1) = tsn(:,:,1,jp_sal) * z3d(:,:,1) 154 CALL iom_put( "sss2" , z3d(:,:,1) ) ! sea surface content of squared salinity 173 155 ELSE 174 CALL iom_put( "toce" , tsn(:,:,:,jp_tem) ) ! temperature 175 IF ( iom_use("sst") ) THEN 176 DO jj = 1, jpj 177 DO ji = 1, jpi 178 z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) 179 END DO 180 END DO 181 CALL iom_put( "sst" , z2d(:,:) ) ! sea surface temperature 182 ENDIF 183 IF ( iom_use("sst2") ) CALL iom_put( "sst2" , z2d(:,:) * z2d(:,:) ) ! square of sea surface temperature 156 CALL iom_put( "toce" , tsn(:,:,:,jp_tem) ) ! temperature 157 CALL iom_put( "sst" , tsn(:,:,1,jp_tem) ) ! sea surface temperature 158 CALL iom_put( "sst2" , tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) ) ! square of sea surface temperature 184 159 CALL iom_put( "soce" , tsn(:,:,:,jp_sal) ) ! salinity 185 IF ( iom_use("sss") ) THEN 186 DO jj = 1, jpj 187 DO ji = 1, jpi 188 z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) 189 END DO 190 END DO 191 CALL iom_put( "sss" , z2d(:,:) ) ! sea surface salinity 192 ENDIF 193 CALL iom_put( "sss2" , z2d(:,:) * z2d(:,:) ) ! square of sea surface salinity 160 CALL iom_put( "sss" , tsn(:,:,1,jp_sal) ) ! sea surface salinity 161 CALL iom_put( "sss2" , tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) ) ! square of sea surface salinity 194 162 END IF 195 163 IF( lk_vvl .AND. (.NOT. ln_dynadv_vec) ) THEN … … 199 167 CALL iom_put( "uoce" , umask(:,:,:) * un(:,:,:) ) ! i-current 200 168 CALL iom_put( "voce" , vmask(:,:,:) * vn(:,:,:) ) ! j-current 201 IF ( iom_use("ssu") ) THEN202 DO jj = 1, jpj203 DO ji = 1, jpi204 z2d(ji,jj) = un(ji,jj,miku(ji,jj))205 END DO206 END DO207 CALL iom_put( "ssu" , z2d ) ! i-current208 ENDIF209 IF ( iom_use("ssv") ) THEN210 DO jj = 1, jpj211 DO ji = 1, jpi212 z2d(ji,jj) = vn(ji,jj,mikv(ji,jj))213 END DO214 END DO215 CALL iom_put( "ssv" , z2d ) ! j-current216 ENDIF217 169 ENDIF 218 170 CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef. -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r4990 r5098 262 262 263 263 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask !: land/ocean mask at T-, U-, V- and F-pts 264 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask !: land/ocean mask at WT-, WU- and WV-pts 264 265 265 266 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: tpol, fpol !: north fold mask (jperio= 3 or 4) … … 332 333 INTEGER FUNCTION dom_oce_alloc() 333 334 !!---------------------------------------------------------------------- 334 INTEGER, DIMENSION(1 1) :: ierr335 INTEGER, DIMENSION(12) :: ierr 335 336 !!---------------------------------------------------------------------- 336 337 ierr(:) = 0 … … 400 401 & vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk), STAT=ierr(11) ) 401 402 403 ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) ) 404 402 405 #if defined key_noslip_accurate 403 ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(1 1) )406 ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(12) ) 404 407 #endif 405 408 ! -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r4990 r5098 281 281 CALL lbc_lnk( fmask_i, 'F', 1._wp ) 282 282 283 ! 3. Ocean/land mask at wu-, wv- and w points 284 !---------------------------------------------- 285 wmask (:,:,1) = tmask(:,:,1) ! ???????? 286 wumask(:,:,1) = umask(:,:,1) ! ???????? 287 wvmask(:,:,1) = vmask(:,:,1) ! ???????? 288 DO jk=2,jpk 289 wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1) 290 wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1) 291 wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1) 292 END DO 283 293 284 294 ! 4. ocean/land mask for the elliptic equation -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r4998 r5098 8 8 !! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: 9 9 !! vvl option includes z_star and z_tilde coordinates 10 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 10 11 !!---------------------------------------------------------------------- 11 12 !! 'key_vvl' variable volume … … 125 126 INTEGER :: ji,jj,jk 126 127 INTEGER :: ii0, ii1, ij0, ij1 128 REAL(wp):: zcoef 127 129 !!---------------------------------------------------------------------- 128 130 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') … … 164 166 ! t- and w- points depth 165 167 ! ---------------------- 168 ! set the isf depth as it is in the initial step 166 169 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 167 170 fsdepw_n(:,:,1) = 0.0_wp … … 169 172 fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1) 170 173 fsdepw_b(:,:,1) = 0.0_wp 171 DO jj = 1,jpj 172 DO ji = 1,jpi 173 DO jk = 2,mikt(ji,jj)-1 174 fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 175 fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 176 fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 177 fsdept_b(ji,jj,jk) = gdept_0(ji,jj,jk) 178 fsdepw_b(ji,jj,jk) = gdepw_0(ji,jj,jk) 179 END DO 180 IF (mikt(ji,jj) .GT. 1) THEN 181 jk = mikt(ji,jj) 182 fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk) 183 fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 184 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk ) - sshn (ji,jj) 185 fsdept_b(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_b(ji,jj,jk) 186 fsdepw_b(ji,jj,jk) = gdepw_0(ji,jj,jk) 187 END IF 188 DO jk = mikt(ji,jj)+1, jpk 189 fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk) 174 175 DO jk = 2, jpk 176 DO jj = 1,jpj 177 DO ji = 1,jpi 178 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 179 ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 180 ! 0.5 where jk = mikt 181 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 190 182 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 191 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk ) - sshn (ji,jj) 192 fsdept_b(ji,jj,jk) = fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk) 183 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 184 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 185 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 193 186 fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1) 187 fsdept_b(ji,jj,jk) = zcoef * ( fsdepw_b(ji,jj,jk ) + 0.5 * fse3w_b(ji,jj,jk)) & 188 & + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk)) 194 189 END DO 195 190 END DO … … 590 585 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def 591 586 INTEGER :: ji,jj,jk ! dummy loop indices 587 REAL(wp) :: zcoef 592 588 !!---------------------------------------------------------------------- 593 589 … … 638 634 ! t- and w- points depth 639 635 ! ---------------------- 636 ! set the isf depth as it is in the initial step 640 637 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 641 638 fsdepw_n(:,:,1) = 0.0_wp 642 639 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 643 DO jj = 1,jpj 644 DO ji = 1,jpi 645 DO jk = 2,mikt(ji,jj)-1 646 fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 647 fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 648 fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 649 END DO 650 IF (mikt(ji,jj) .GT. 1) THEN 651 jk = mikt(ji,jj) 652 fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk) 653 fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 654 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk ) - sshn (ji,jj) 655 END IF 656 DO jk = mikt(ji,jj)+1, jpk 657 fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk) 640 641 DO jk = 2, jpk 642 DO jj = 1,jpj 643 DO ji = 1,jpi 644 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 645 ! 1 for jk = mikt 646 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 658 647 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 659 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk ) - sshn (ji,jj) 648 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 649 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 650 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 660 651 END DO 661 652 END DO 662 653 END DO 654 663 655 ! Local depth and Inverse of the local depth of the water column at u- and v- points 664 656 ! ---------------------------------------------------------------------------------- -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5040 r5098 17 17 !! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function 18 18 !! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case 19 !! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye 19 20 !!---------------------------------------------------------------------- 20 21 … … 35 36 USE oce ! ocean variables 36 37 USE dom_oce ! ocean domain 37 USE sbc_oce ! surface variable (isf)38 38 USE closea ! closed seas 39 39 USE c1d ! 1D vertical configuration … … 472 472 ! 473 473 ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 474 IF( cp_cfg == "isomip" ) THEN474 IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN 475 475 ! 476 476 risfdep(:,:)=200.e0 … … 482 482 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 483 483 ! 484 ELSEIF ( cp_cfg == "isomip2" ) THEN484 ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN 485 485 ! 486 486 risfdep(:,:)=0.e0 … … 584 584 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 585 585 ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin 586 WHERE (bathy == risfdep) 587 bathy = 0.0_wp ; risfdep = 0.0_wp 588 END WHERE 586 IF ( ln_isfcav ) THEN 587 WHERE (bathy == risfdep) 588 bathy = 0.0_wp ; risfdep = 0.0_wp 589 END WHERE 590 END IF 589 591 ! end patch 590 592 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level … … 964 966 INTEGER :: ik, it ! temporary integers 965 967 INTEGER :: id, jd, nprocd 968 LOGICAL :: ll_print ! Allow control print for debugging 969 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 970 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 971 REAL(wp) :: zmax, zmin ! Maximum and minimum depth 972 REAL(wp) :: zdiff ! temporary scalar 973 REAL(wp) :: zrefdep ! temporary scalar 974 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 975 !!--------------------------------------------------------------------- 976 ! 977 IF( nn_timing == 1 ) CALL timing_start('zgr_zps') 978 ! 979 CALL wrk_alloc( jpi, jpj, jpk, zprt ) 980 ! 981 IF(lwp) WRITE(numout,*) 982 IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps' 983 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 984 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 985 986 ll_print = .FALSE. ! Local variable for debugging 987 988 IF(lwp .AND. ll_print) THEN ! control print of the ocean depth 989 WRITE(numout,*) 990 WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)' 991 CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 992 ENDIF 993 994 995 ! bathymetry in level (from bathy_meter) 996 ! =================== 997 zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 998 bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat) 999 WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0 1000 ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level 1001 END WHERE 1002 1003 ! Compute mbathy for ocean points (i.e. the number of ocean levels) 1004 ! find the number of ocean levels such that the last level thickness 1005 ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 1006 ! e3t_1d is the reference level thickness 1007 DO jk = jpkm1, 1, -1 1008 zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1009 WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 1010 END DO 1011 1012 IF ( ln_isfcav ) CALL zgr_isf 1013 1014 ! Scale factors and depth at T- and W-points 1015 DO jk = 1, jpk ! intitialization to the reference z-coordinate 1016 gdept_0(:,:,jk) = gdept_1d(jk) 1017 gdepw_0(:,:,jk) = gdepw_1d(jk) 1018 e3t_0 (:,:,jk) = e3t_1d (jk) 1019 e3w_0 (:,:,jk) = e3w_1d (jk) 1020 END DO 1021 ! 1022 DO jj = 1, jpj 1023 DO ji = 1, jpi 1024 ik = mbathy(ji,jj) 1025 IF( ik > 0 ) THEN ! ocean point only 1026 ! max ocean level case 1027 IF( ik == jpkm1 ) THEN 1028 zdepwp = bathy(ji,jj) 1029 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1030 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1031 e3t_0(ji,jj,ik ) = ze3tp 1032 e3t_0(ji,jj,ik+1) = ze3tp 1033 e3w_0(ji,jj,ik ) = ze3wp 1034 e3w_0(ji,jj,ik+1) = ze3tp 1035 gdepw_0(ji,jj,ik+1) = zdepwp 1036 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1037 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1038 ! 1039 ELSE ! standard case 1040 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1041 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1042 ENDIF 1043 !gm Bug? check the gdepw_1d 1044 ! ... on ik 1045 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1046 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1047 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1048 e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1049 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1050 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1051 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1052 ! ... on ik+1 1053 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1054 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1055 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 1056 ENDIF 1057 ENDIF 1058 END DO 1059 END DO 1060 ! 1061 it = 0 1062 DO jj = 1, jpj 1063 DO ji = 1, jpi 1064 ik = mbathy(ji,jj) 1065 IF( ik > 0 ) THEN ! ocean point only 1066 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1067 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1068 ! test 1069 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1070 IF( zdiff <= 0._wp .AND. lwp ) THEN 1071 it = it + 1 1072 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1073 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1074 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1075 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1076 ENDIF 1077 ENDIF 1078 END DO 1079 END DO 1080 ! 1081 IF ( ln_isfcav ) THEN 1082 ! (ISF) Definition of e3t, u, v, w for ISF case 1083 DO jj = 1, jpj 1084 DO ji = 1, jpi 1085 ik = misfdep(ji,jj) 1086 IF( ik > 1 ) THEN ! ice shelf point only 1087 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1088 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1089 !gm Bug? check the gdepw_0 1090 ! ... on ik 1091 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1092 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1093 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1094 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1095 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1096 1097 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1098 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1099 ENDIF 1100 ! ... on ik / ik-1 1101 e3w_0 (ji,jj,ik ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1102 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1103 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1104 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1105 ENDIF 1106 END DO 1107 END DO 1108 ! 1109 it = 0 1110 DO jj = 1, jpj 1111 DO ji = 1, jpi 1112 ik = misfdep(ji,jj) 1113 IF( ik > 1 ) THEN ! ice shelf point only 1114 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1115 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1116 ! test 1117 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1118 IF( zdiff <= 0. .AND. lwp ) THEN 1119 it = it + 1 1120 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1121 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1122 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1123 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1124 ENDIF 1125 ENDIF 1126 END DO 1127 END DO 1128 END IF 1129 ! END (ISF) 1130 1131 ! Scale factors and depth at U-, V-, UW and VW-points 1132 DO jk = 1, jpk ! initialisation to z-scale factors 1133 e3u_0 (:,:,jk) = e3t_1d(jk) 1134 e3v_0 (:,:,jk) = e3t_1d(jk) 1135 e3uw_0(:,:,jk) = e3w_1d(jk) 1136 e3vw_0(:,:,jk) = e3w_1d(jk) 1137 END DO 1138 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 1139 DO jj = 1, jpjm1 1140 DO ji = 1, fs_jpim1 ! vector opt. 1141 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 1142 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 1143 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 1144 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 1145 END DO 1146 END DO 1147 END DO 1148 !IF ( ln_isfcav ) THEN 1149 ! (ISF) define e3uw (adapted for 2 cells in the water column) 1150 DO jk = 2,jpk 1151 DO jj = 1, jpjm1 1152 DO ji = 1, fs_jpim1 ! vector opt. 1153 e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj ,jk) ) & 1154 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj ,jk-1) ) 1155 e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji ,jj+1,jk) ) & 1156 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji ,jj+1,jk-1) ) 1157 END DO 1158 END DO 1159 END DO 1160 !END IF 1161 !End (ISF) 1162 1163 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 1164 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 1165 ! 1166 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1167 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk) 1168 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk) 1169 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk) 1170 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk) 1171 END DO 1172 1173 ! Scale factor at F-point 1174 DO jk = 1, jpk ! initialisation to z-scale factors 1175 e3f_0(:,:,jk) = e3t_1d(jk) 1176 END DO 1177 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors 1178 DO jj = 1, jpjm1 1179 DO ji = 1, fs_jpim1 ! vector opt. 1180 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 1181 END DO 1182 END DO 1183 END DO 1184 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions 1185 ! 1186 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1187 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk) 1188 END DO 1189 !!gm bug ? : must be a do loop with mj0,mj1 1190 ! 1191 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 1192 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 1193 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 1194 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 1195 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 1196 1197 ! Control of the sign 1198 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' ) 1199 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' ) 1200 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' ) 1201 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' ) 1202 1203 ! Compute gdep3w_0 (vertical sum of e3w) 1204 IF ( ln_isfcav ) THEN 1205 WHERE (misfdep == 0) misfdep = 1 1206 DO jj = 1,jpj 1207 DO ji = 1,jpi 1208 gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1209 DO jk = 2, misfdep(ji,jj) 1210 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1211 END DO 1212 IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 1213 DO jk = misfdep(ji,jj) + 1, jpk 1214 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1215 END DO 1216 END DO 1217 END DO 1218 ELSE ! no cavity 1219 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1220 DO jk = 2, jpk 1221 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1222 END DO 1223 END IF 1224 ! ! ================= ! 1225 IF(lwp .AND. ll_print) THEN ! Control print ! 1226 ! ! ================= ! 1227 DO jj = 1,jpj 1228 DO ji = 1, jpi 1229 ik = MAX( mbathy(ji,jj), 1 ) 1230 zprt(ji,jj,1) = e3t_0 (ji,jj,ik) 1231 zprt(ji,jj,2) = e3w_0 (ji,jj,ik) 1232 zprt(ji,jj,3) = e3u_0 (ji,jj,ik) 1233 zprt(ji,jj,4) = e3v_0 (ji,jj,ik) 1234 zprt(ji,jj,5) = e3f_0 (ji,jj,ik) 1235 zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 1236 END DO 1237 END DO 1238 WRITE(numout,*) 1239 WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1240 WRITE(numout,*) 1241 WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1242 WRITE(numout,*) 1243 WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1244 WRITE(numout,*) 1245 WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1246 WRITE(numout,*) 1247 WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1248 WRITE(numout,*) 1249 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1250 ENDIF 1251 ! 1252 CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 1253 ! 1254 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps') 1255 ! 1256 END SUBROUTINE zgr_zps 1257 1258 SUBROUTINE zgr_isf 1259 !!---------------------------------------------------------------------- 1260 !! *** ROUTINE zgr_bat_ctl *** 1261 !! 1262 !! ** Purpose : check the bathymetry in levels 1263 !! 1264 !! ** Method : THe water column have to contained at least 2 cells 1265 !! Bathymetry and isfdraft are modified (dig/close) to respect 1266 !! this criterion. 1267 !! 1268 !! 1269 !! ** Action : - test compatibility between isfdraft and bathy 1270 !! - bathy and isfdraft are modified 1271 !!---------------------------------------------------------------------- 1272 !! 1273 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1274 INTEGER :: ik, it ! temporary integers 1275 INTEGER :: id, jd, nprocd 966 1276 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF) 967 1277 LOGICAL :: ll_print ! Allow control print for debugging … … 971 1281 REAL(wp) :: zdiff ! temporary scalar 972 1282 REAL(wp) :: zrefdep ! temporary scalar 973 REAL(wp) :: zbathydiff, zrisfdepdiff 974 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 3D workspace (ISH) 975 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 3D workspace (ISH) 976 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 1283 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar 1284 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH) 1285 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH) 977 1286 !!--------------------------------------------------------------------- 978 1287 ! 979 IF( nn_timing == 1 ) CALL timing_start('zgr_zps') 980 ! 981 CALL wrk_alloc( jpi, jpj, jpk, zprt ) 1288 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1289 ! 982 1290 CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 983 CALL wrk_alloc( jpi, jpj, zmbathy, zmisfdep) 984 ! 985 IF(lwp) WRITE(numout,*) 986 IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps' 987 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 988 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 989 990 ll_print = .FALSE. ! Local variable for debugging 991 992 IF(lwp .AND. ll_print) THEN ! control print of the ocean depth 993 WRITE(numout,*) 994 WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)' 995 CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 996 ENDIF 997 998 ! bathymetry in level (from bathy_meter) 999 ! =================== 1000 zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 1001 bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat) 1002 WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0 1003 ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level 1004 END WHERE 1005 1006 ! Compute mbathy for ocean points (i.e. the number of ocean levels) 1007 ! find the number of ocean levels such that the last level thickness 1008 ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 1009 ! e3t_1d is the reference level thickness 1010 DO jk = jpkm1, 1, -1 1011 zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1012 WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 1013 END DO 1291 CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy ) 1292 1293 1014 1294 ! (ISF) compute misfdep 1015 1295 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ; misfdep(:,:) = 1 ! open water : set misfdep to 1 … … 1475 1755 ENDIF 1476 1756 1477 ! Scale factors and depth at T- and W-points1478 DO jk = 1, jpk ! intitialization to the reference z-coordinate1479 gdept_0(:,:,jk) = gdept_1d(jk)1480 gdepw_0(:,:,jk) = gdepw_1d(jk)1481 e3t_0 (:,:,jk) = e3t_1d (jk)1482 e3w_0 (:,:,jk) = e3w_1d (jk)1483 END DO1484 !1485 DO jj = 1, jpj1486 DO ji = 1, jpi1487 ik = mbathy(ji,jj)1488 IF( ik > 0 ) THEN ! ocean point only1489 ! max ocean level case1490 IF( ik == jpkm1 ) THEN1491 zdepwp = bathy(ji,jj)1492 ze3tp = bathy(ji,jj) - gdepw_1d(ik)1493 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )1494 e3t_0(ji,jj,ik ) = ze3tp1495 e3t_0(ji,jj,ik+1) = ze3tp1496 e3w_0(ji,jj,ik ) = ze3wp1497 e3w_0(ji,jj,ik+1) = ze3tp1498 gdepw_0(ji,jj,ik+1) = zdepwp1499 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp1500 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp1501 !1502 ELSE ! standard case1503 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)1504 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)1505 ENDIF1506 !gm Bug? check the gdepw_1d1507 ! ... on ik1508 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &1509 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) &1510 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))1511 e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &1512 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )1513 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &1514 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )1515 ! ... on ik+11516 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)1517 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)1518 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)1519 ENDIF1520 ENDIF1521 END DO1522 END DO1523 !1524 it = 01525 DO jj = 1, jpj1526 DO ji = 1, jpi1527 ik = mbathy(ji,jj)1528 IF( ik > 0 ) THEN ! ocean point only1529 e3tp (ji,jj) = e3t_0(ji,jj,ik)1530 e3wp (ji,jj) = e3w_0(ji,jj,ik)1531 ! test1532 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )1533 IF( zdiff <= 0._wp .AND. lwp ) THEN1534 it = it + 11535 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj1536 WRITE(numout,*) ' bathy = ', bathy(ji,jj)1537 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff1538 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )1539 ENDIF1540 ENDIF1541 END DO1542 END DO1543 !1544 ! (ISF) Definition of e3t, u, v, w for ISF case1545 DO jj = 1, jpj1546 DO ji = 1, jpi1547 ik = misfdep(ji,jj)1548 IF( ik > 1 ) THEN ! ice shelf point only1549 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik)1550 gdepw_0(ji,jj,ik) = risfdep(ji,jj)1551 !gm Bug? check the gdepw_01552 ! ... on ik1553 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) &1554 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) &1555 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) )1556 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)1557 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)1558 1559 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column)1560 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)1561 ENDIF1562 ! ... on ik / ik-11563 e3w_0 (ji,jj,ik ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))1564 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)1565 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code1566 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)1567 ENDIF1568 END DO1569 END DO1570 !1571 it = 01572 DO jj = 1, jpj1573 DO ji = 1, jpi1574 ik = misfdep(ji,jj)1575 IF( ik > 1 ) THEN ! ice shelf point only1576 e3tp (ji,jj) = e3t_0(ji,jj,ik )1577 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )1578 ! test1579 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik )1580 IF( zdiff <= 0. .AND. lwp ) THEN1581 it = it + 11582 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj1583 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)1584 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff1585 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj)1586 ENDIF1587 ENDIF1588 END DO1589 END DO1590 ! END (ISF)1591 1592 ! Scale factors and depth at U-, V-, UW and VW-points1593 DO jk = 1, jpk ! initialisation to z-scale factors1594 e3u_0 (:,:,jk) = e3t_1d(jk)1595 e3v_0 (:,:,jk) = e3t_1d(jk)1596 e3uw_0(:,:,jk) = e3w_1d(jk)1597 e3vw_0(:,:,jk) = e3w_1d(jk)1598 END DO1599 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors1600 DO jj = 1, jpjm11601 DO ji = 1, fs_jpim1 ! vector opt.1602 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )1603 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )1604 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )1605 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )1606 END DO1607 END DO1608 END DO1609 ! (ISF) define e3uw1610 DO jk = 2,jpk1611 DO jj = 1, jpjm11612 DO ji = 1, fs_jpim1 ! vector opt.1613 e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj ,jk) ) &1614 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj ,jk-1) )1615 e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji ,jj+1,jk) ) &1616 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji ,jj+1,jk-1) )1617 END DO1618 END DO1619 END DO1620 !End (ISF)1621 1622 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions1623 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp )1624 !1625 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1626 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)1627 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)1628 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)1629 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)1630 END DO1631 1632 ! Scale factor at F-point1633 DO jk = 1, jpk ! initialisation to z-scale factors1634 e3f_0(:,:,jk) = e3t_1d(jk)1635 END DO1636 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors1637 DO jj = 1, jpjm11638 DO ji = 1, fs_jpim1 ! vector opt.1639 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )1640 END DO1641 END DO1642 END DO1643 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions1644 !1645 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1646 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)1647 END DO1648 !!gm bug ? : must be a do loop with mj0,mj11649 !1650 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 21651 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)1652 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)1653 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)1654 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)1655 1656 ! Control of the sign1657 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )1658 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )1659 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )1660 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )1661 1662 ! Compute gdep3w_0 (vertical sum of e3w)1663 WHERE (misfdep == 0) misfdep = 11664 DO jj = 1,jpj1665 DO ji = 1,jpi1666 gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)1667 DO jk = 2, misfdep(ji,jj)1668 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1669 END DO1670 IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))1671 DO jk = misfdep(ji,jj) + 1, jpk1672 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1673 END DO1674 END DO1675 END DO1676 ! ! ================= !1677 IF(lwp .AND. ll_print) THEN ! Control print !1678 ! ! ================= !1679 DO jj = 1,jpj1680 DO ji = 1, jpi1681 ik = MAX( mbathy(ji,jj), 1 )1682 zprt(ji,jj,1) = e3t_0 (ji,jj,ik)1683 zprt(ji,jj,2) = e3w_0 (ji,jj,ik)1684 zprt(ji,jj,3) = e3u_0 (ji,jj,ik)1685 zprt(ji,jj,4) = e3v_0 (ji,jj,ik)1686 zprt(ji,jj,5) = e3f_0 (ji,jj,ik)1687 zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)1688 END DO1689 END DO1690 WRITE(numout,*)1691 WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1692 WRITE(numout,*)1693 WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1694 WRITE(numout,*)1695 WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1696 WRITE(numout,*)1697 WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1698 WRITE(numout,*)1699 WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1700 WRITE(numout,*)1701 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1702 ENDIF1703 !1704 CALL wrk_dealloc( jpi, jpj, jpk, zprt )1705 1757 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1706 1758 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1707 ! 1708 IF( nn_timing == 1 ) CALL timing_stop('zgr_ zps')1709 !1710 END SUBROUTINE zgr_zps1759 1760 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1761 1762 END SUBROUTINE 1711 1763 1712 1764 SUBROUTINE zgr_sco -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r4990 r5098 137 137 CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! before potential and in situ densities 138 138 #if ! defined key_c1d 139 IF( ln_zps ) CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 140 & rhd, gru , grv, aru, arv, gzu, gzv, ge3ru, ge3rv, & ! 141 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 139 IF( ln_zps ) CALL zps_hde ( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 140 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 141 IF( ln_zps .AND. ln_isfcav) & 142 & CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps for top cell (ISF) 143 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 144 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 142 145 #endif 143 146 ! -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90
r4990 r5098 17 17 !! 3.3 ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module 18 18 !! - ! 2010-10 (R. Furner, G. Madec) runoff and cla added directly here 19 !! 3.6 ! 2014-11 (P. Mathiot) isf added directly here 19 20 !!---------------------------------------------------------------------- 20 21 -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
r4990 r5098 80 80 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu) 81 81 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv) 82 83 ! (ISF) stability criteria for top friction84 ikbu = miku(ji,jj) ! first wet ocean u- & v-levels85 ikbv = mikv(ji,jj)86 !87 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)88 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( tfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu) &89 & * (1.-umask(ji,jj,1))90 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( tfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv) &91 & * (1.-vmask(ji,jj,1))92 ! (ISF)93 94 82 END DO 95 83 END DO 84 85 IF ( ln_isfcav ) THEN 86 DO jj = 2, jpjm1 87 DO ji = 2, jpim1 88 ! (ISF) stability criteria for top friction 89 ikbu = miku(ji,jj) ! first wet ocean u- & v-levels 90 ikbv = mikv(ji,jj) 91 ! 92 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 93 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( tfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu) & 94 & * (1.-umask(ji,jj,1)) 95 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( tfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv) & 96 & * (1.-vmask(ji,jj,1)) 97 ! (ISF) 98 END DO 99 END DO 100 END IF 96 101 97 102 ! -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r4990 r5098 16 16 !! 3.4 ! 2011-11 (H. Liu) hpg_prj: Original code for s-coordinates 17 17 !! ! (A. Coward) suppression of hel, wdj and rot options 18 !! 3.6 ! 2014-11 (P. Mathiot) hpg_isf: original code for ice shelf cavity 18 19 !!---------------------------------------------------------------------- 19 20 … … 25 26 !! hpg_zps : z-coordinate plus partial steps (interpolation) 26 27 !! hpg_sco : s-coordinate (standard jacobian formulation) 28 !! hpg_isf : s-coordinate (sco formulation) adapted to ice shelf 27 29 !! hpg_djc : s-coordinate (Density Jacobian with Cubic polynomial) 28 30 !! hpg_prj : s-coordinate (Pressure Jacobian with Cubic polynomial) … … 55 57 LOGICAL , PUBLIC :: ln_hpg_djc !: s-coordinate (Density Jacobian with Cubic polynomial) 56 58 LOGICAL , PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) 59 LOGICAL , PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf 57 60 LOGICAL , PUBLIC :: ln_dynhpg_imp !: semi-implicite hpg flag 58 61 … … 97 100 CASE ( 3 ) ; CALL hpg_djc ( kt ) ! s-coordinate (Density Jacobian with Cubic polynomial) 98 101 CASE ( 4 ) ; CALL hpg_prj ( kt ) ! s-coordinate (Pressure Jacobian scheme) 102 CASE ( 5 ) ; CALL hpg_isf ( kt ) ! s-coordinate similar to sco modify for ice shelf 99 103 END SELECT 100 104 ! … … 128 132 !! 129 133 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 130 & ln_hpg_djc, ln_hpg_prj, ln_ dynhpg_imp134 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf, ln_dynhpg_imp 131 135 !!---------------------------------------------------------------------- 132 136 ! … … 148 152 WRITE(numout,*) ' z-coord. - partial steps (interpolation) ln_hpg_zps = ', ln_hpg_zps 149 153 WRITE(numout,*) ' s-coord. (standard jacobian formulation) ln_hpg_sco = ', ln_hpg_sco 154 WRITE(numout,*) ' s-coord. (standard jacobian formulation) for isf ln_hpg_isf = ', ln_hpg_isf 150 155 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 151 156 WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj … … 158 163 & either ln_hpg_sco or ln_hpg_prj instead') 159 164 ! 160 IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj ) ) &165 IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) ) & 161 166 & CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:& 162 167 & the standard jacobian formulation hpg_sco or & … … 169 174 IF( ln_hpg_djc ) nhpg = 3 170 175 IF( ln_hpg_prj ) nhpg = 4 176 IF( ln_hpg_isf ) nhpg = 5 171 177 ! 172 178 ! ! Consistency check … … 177 183 IF( ln_hpg_djc ) ioptio = ioptio + 1 178 184 IF( ln_hpg_prj ) ioptio = ioptio + 1 185 IF( ln_hpg_isf ) ioptio = ioptio + 1 179 186 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 180 IF( ( ln_hpg_zco .OR. ln_hpg_zps .OR. ln_hpg_djc .OR. ln_hpg_prj ) .AND. nn_isf .NE. 0) &181 & CALL ctl_stop( 'Only hpg_ scohas been corrected to work with ice shelf cavity.' )187 IF( ( .NOT. ln_hpg_isf ) .AND. ln_isfcav ) & 188 & CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' ) 182 189 ! 183 190 END SUBROUTINE dyn_hpg_init … … 345 352 END SUBROUTINE hpg_zps 346 353 347 348 354 SUBROUTINE hpg_sco( kt ) 349 355 !!--------------------------------------------------------------------- … … 366 372 INTEGER, INTENT(in) :: kt ! ocean time-step index 367 373 !! 374 INTEGER :: ji, jj, jk ! dummy loop indices 375 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 376 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 377 !!---------------------------------------------------------------------- 378 ! 379 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 380 ! 381 IF( kt == nit000 ) THEN 382 IF(lwp) WRITE(numout,*) 383 IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' 384 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 385 ENDIF 386 387 ! Local constant initialization 388 zcoef0 = - grav * 0.5_wp 389 ! To use density and not density anomaly 390 IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume 391 ELSE ; znad = 0._wp ! Fixed volume 392 ENDIF 393 394 ! Surface value 395 DO jj = 2, jpjm1 396 DO ji = fs_2, fs_jpim1 ! vector opt. 397 ! hydrostatic pressure gradient along s-surfaces 398 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & 399 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 400 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & 401 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 402 ! s-coordinate pressure gradient correction 403 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 404 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 405 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 406 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 407 ! add to the general momentum trend 408 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 409 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 410 END DO 411 END DO 412 413 ! interior value (2=<jk=<jpkm1) 414 DO jk = 2, jpkm1 415 DO jj = 2, jpjm1 416 DO ji = fs_2, fs_jpim1 ! vector opt. 417 ! hydrostatic pressure gradient along s-surfaces 418 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 419 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 420 & - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 421 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 422 & * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 423 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 424 ! s-coordinate pressure gradient correction 425 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 426 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 427 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 428 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 429 ! add to the general momentum trend 430 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 431 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 432 END DO 433 END DO 434 END DO 435 ! 436 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 437 ! 438 END SUBROUTINE hpg_sco 439 440 SUBROUTINE hpg_isf( kt ) 441 !!--------------------------------------------------------------------- 442 !! *** ROUTINE hpg_sco *** 443 !! 444 !! ** Method : s-coordinate case. Jacobian scheme. 445 !! The now hydrostatic pressure gradient at a given level, jk, 446 !! is computed by taking the vertical integral of the in-situ 447 !! density gradient along the model level from the suface to that 448 !! level. s-coordinates (ln_sco): a corrective term is added 449 !! to the horizontal pressure gradient : 450 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] 451 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ] 452 !! add it to the general momentum trend (ua,va). 453 !! ua = ua - 1/e1u * zhpi 454 !! va = va - 1/e2v * zhpj 455 !! iceload is added and partial cell case are added to the top and bottom 456 !! 457 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 458 !!---------------------------------------------------------------------- 459 INTEGER, INTENT(in) :: kt ! ocean time-step index 460 !! 368 461 INTEGER :: ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j ! dummy loop indices 369 462 REAL(wp) :: zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1 ! temporary scalars … … 379 472 IF( kt == nit000 ) THEN 380 473 IF(lwp) WRITE(numout,*) 381 IF(lwp) WRITE(numout,*) 'dyn:hpg_ sco : hydrostatic pressure gradient trend'474 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 382 475 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 383 476 ENDIF … … 565 658 !================================================================================== 566 659 567 # if defined key_vectopt_loop568 jj = 1569 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)570 # else571 660 DO jj = 2, jpjm1 572 661 DO ji = 2, jpim1 573 # endif574 662 iku = mbku(ji,jj) 575 663 ikv = mbkv(ji,jj) … … 598 686 va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 599 687 END IF 600 # if ! defined key_vectopt_loop 601 END DO 602 # endif 688 END DO 603 689 END DO 604 690 … … 610 696 CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 611 697 ! 612 END SUBROUTINE hpg_ sco698 END SUBROUTINE hpg_isf 613 699 614 700 -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r4990 r5098 250 250 IF( ( ioptio > 1 .AND. .NOT. lk_esopa ) .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) ) & 251 251 & CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' ) 252 IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. nn_isf .NE. 0) &252 IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. ln_isfcav ) & 253 253 & CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' ) 254 254 ! -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5032 r5098 22 22 USE dom_oce ! ocean space and time domain 23 23 USE sbc_oce ! surface boundary condition: ocean 24 USE sbcisf ! ice shelf variable (fwfisf) 24 25 USE dynspg_oce ! surface pressure gradient variables 25 26 USE phycst ! physical constants … … 453 454 ! ! Surface net water flux and rivers 454 455 IF (ln_bt_fw) THEN 455 zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) )456 zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) ) 456 457 ELSE 457 zssh_frc(:,:) = zraur * z1_2 * (emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)) 458 zssh_frc(:,:) = zraur * z1_2 * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) & 459 & + rdivisf * ( fwfisf(:,:) + fwfisf_b(:,:) ) ) 458 460 ENDIF 459 461 #if defined key_asminc -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
r4990 r5098 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) ) !* wumask(ji,jj,jk) 93 zwvw(ji,jj,jk) = ( zww(ji ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1)-vn(ji,jj,jk) ) !* wvmask(ji,jj,jk) 94 94 END DO 95 95 END DO … … 228 228 DO jj = 2, jpjm1 ! vertical momentum advection at w-point 229 229 DO ji = fs_2, fs_jpim1 ! vector opt. 230 zwuw(ji,jj,jk) = ( zww(ji+1,jj ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) 231 zwvw(ji,jj,jk) = ( zww(ji ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) 230 zwuw(ji,jj,jk) = ( zww(ji+1,jj ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) !* wumask(ji,jj,jk) 231 zwvw(ji,jj,jk) = ( zww(ji ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) !* wvmask(ji,jj,jk) 232 232 END DO 233 233 END DO -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r4990 r5098 105 105 avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 106 106 avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 107 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 108 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 109 IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu) 110 IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv) 111 END DO 112 END DO 107 END DO 108 END DO 109 IF ( ln_isfcav ) THEN 110 DO jj = 2, jpjm1 111 DO ji = 2, jpim1 112 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 113 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 114 IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu) 115 IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv) 116 END DO 117 END DO 118 END IF 113 119 ENDIF 114 120 … … 168 174 zcoef = - p2dt / ze3ua 169 175 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk ) 170 zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk)176 zwi(ji,jj,jk) = zzwi * wumask(ji,jj,jk) 171 177 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 172 zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1)173 zwd(ji,jj,jk) = 1._wp - z wi(ji,jj,jk)- zzws178 zws(ji,jj,jk) = zzws * wumask(ji,jj,jk+1) 179 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 174 180 END DO 175 181 END DO … … 198 204 ! 199 205 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 200 DO j j = 2, jpjm1201 DO j i = fs_2, fs_jpim1 ! vector opt.202 DO j k = miku(ji,jj)+1, jpkm1206 DO jk = 2, jpkm1 207 DO jj = 2, jpjm1 208 DO ji = fs_2, fs_jpim1 ! vector opt. 203 209 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 204 210 END DO … … 208 214 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 209 215 DO ji = fs_2, fs_jpim1 ! vector opt. 210 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,miku(ji,jj)) + r_vvl * fse3u_a(ji,jj,miku(ji,jj))211 216 #if defined key_dynspg_ts 212 ua(ji,jj,miku(ji,jj)) = ua(ji,jj,miku(ji,jj)) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 213 & / ( ze3ua * rau0 ) 217 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl * fse3u_a(ji,jj,1) 218 ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 219 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 214 220 #else 215 ua(ji,jj,miku(ji,jj)) = ub(ji,jj,miku(ji,jj)) & 216 & + p2dt *(ua(ji,jj,miku(ji,jj)) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 217 & / ( fse3u(ji,jj,miku(ji,jj)) * rau0 ) ) 218 #endif 219 DO jk = miku(ji,jj)+1, jpkm1 221 ua(ji,jj,1) = ub(ji,jj,1) & 222 & + p2dt *(ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 223 & / ( fse3u(ji,jj,1) * rau0 ) * umask(ji,jj,1) ) 224 #endif 225 END DO 226 END DO 227 DO jk = 2, jpkm1 228 DO jj = 2, jpjm1 229 DO ji = fs_2, fs_jpim1 220 230 #if defined key_dynspg_ts 221 231 zrhs = ua(ji,jj,jk) ! zrhs=right hand side … … 231 241 DO ji = fs_2, fs_jpim1 ! vector opt. 232 242 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 233 DO jk = jpk-2, miku(ji,jj), -1 243 END DO 244 END DO 245 DO jk = jpk-2, 1, -1 246 DO jj = 2, jpjm1 247 DO ji = fs_2, fs_jpim1 234 248 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 235 249 END DO … … 260 274 zcoef = - p2dt / ze3va 261 275 zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk ) 262 zwi(ji,jj,jk) = zzwi * vmask(ji,jj,jk)276 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk) 263 277 zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 264 zws(ji,jj,jk) = zzws * vmask(ji,jj,jk+1)265 zwd(ji,jj,jk) = 1._wp - z wi(ji,jj,jk)- zzws278 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 279 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 266 280 END DO 267 281 END DO … … 290 304 ! 291 305 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 292 DO j j = 2, jpjm1293 DO j i = fs_2, fs_jpim1 ! vector opt.294 DO j k = mikv(ji,jj)+1, jpkm1306 DO jk = 2, jpkm1 307 DO jj = 2, jpjm1 308 DO ji = fs_2, fs_jpim1 ! vector opt. 295 309 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 296 310 END DO … … 300 314 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 301 315 DO ji = fs_2, fs_jpim1 ! vector opt. 302 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,mikv(ji,jj)) + r_vvl * fse3v_a(ji,jj,mikv(ji,jj))303 316 #if defined key_dynspg_ts 304 va(ji,jj,mikv(ji,jj)) = va(ji,jj,mikv(ji,jj)) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 317 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl * fse3v_a(ji,jj,1) 318 va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 305 319 & / ( ze3va * rau0 ) 306 320 #else 307 va(ji,jj,mikv(ji,jj)) = vb(ji,jj,mikv(ji,jj)) & 308 & + p2dt *(va(ji,jj,mikv(ji,jj)) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 309 & / ( fse3v(ji,jj,mikv(ji,jj)) * rau0 ) ) 310 #endif 311 DO jk = mikv(ji,jj)+1, jpkm1 321 va(ji,jj,1) = vb(ji,jj,1) & 322 & + p2dt *(va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 323 & / ( fse3v(ji,jj,1) * rau0 ) ) 324 #endif 325 END DO 326 END DO 327 DO jk = 2, jpkm1 328 DO jj = 2, jpjm1 329 DO ji = fs_2, fs_jpim1 ! vector opt. 312 330 #if defined key_dynspg_ts 313 331 zrhs = va(ji,jj,jk) ! zrhs=right hand side … … 323 341 DO ji = fs_2, fs_jpim1 ! vector opt. 324 342 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 325 DO jk = jpk-2, mikv(ji,jj), -1 343 END DO 344 END DO 345 DO jk = jpk-2, 1, -1 346 DO jj = 2, jpjm1 347 DO ji = fs_2, fs_jpim1 326 348 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 327 349 END DO … … 349 371 avmu(ji,jj,ikbu+1) = 0.e0 350 372 avmv(ji,jj,ikbv+1) = 0.e0 351 ikbu = miku(ji,jj) ! ocean top level at u- and v-points352 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points)353 IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0354 IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0355 373 END DO 356 374 END DO 375 IF (ln_isfcav) THEN 376 DO jj = 2, jpjm1 377 DO ji = 2, jpim1 378 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 379 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 380 IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0 381 IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0 382 END DO 383 END DO 384 END IF 357 385 ENDIF 358 386 ! -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r5016 r5098 143 143 DO ji = 1, jpim1 144 144 ! IF should be useless check zpshde (PM) 145 IF ( mbku(ji,jj) > 1 ) zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 146 IF ( mbkv(ji,jj) > 1 ) zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 145 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 146 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 147 END DO 148 END DO 149 ENDIF 150 IF( ln_zps .AND. ln_isfcav ) THEN ! partial steps correction at the bottom ocean level 151 DO jj = 1, jpjm1 152 DO ji = 1, jpim1 147 153 IF ( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj) 148 154 IF ( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj) … … 151 157 ENDIF 152 158 ! 153 zdzr(:,:,1) = 0._wp !== Local vertical density gradient at T-point == ! (evaluated from N^2) 154 DO jk = 1, jpkm1 159 !== Local vertical density gradient at T-point == ! (evaluated from N^2) 160 ! surface initialisation 161 zdzr(:,:,1) = 0._wp 162 ! interior value 163 DO jk = 2, jpkm1 155 164 ! ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 156 165 ! ! trick: tmask(ik ) = 0 => all pn2 = 0 => zdzr = 0 … … 161 170 & * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) ) 162 171 END DO 163 ! surface initialisation 164 DO jj = 1, jpjm1 165 DO ji = 1, jpim1 166 zdzr(ji,jj,1:mikt(ji,jj)) = 0._wp 167 END DO 168 END DO 172 IF ( ln_isfcav ) THEN 173 ! if isf need to overwrite the interior value at at the first ocean point 174 DO jj = 1, jpjm1 175 DO ji = 1, jpim1 176 zdzr(ji,jj,1:mikt(ji,jj)) = 0._wp 177 END DO 178 END DO 179 END IF 169 180 ! 170 181 ! !== Slopes just below the mixed layer ==! … … 210 221 & + zfi * uslpml(ji,jj) & 211 222 & * zdepu / MAX( zhmlpu(ji,jj), 5._wp ) 212 zwz(ji,jj,jk) = zwz(ji,jj,jk) * umask(ji,jj,jk) * umask(ji,jj,jk-1)223 zwz(ji,jj,jk) = zwz(ji,jj,jk) * wumask(ji,jj,jk) 213 224 zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps ) & 214 225 & + zfj * vslpml(ji,jj) & 215 226 & * zdepv / MAX( zhmlpv(ji,jj), 5._wp ) 216 zww(ji,jj,jk) = zww(ji,jj,jk) * vmask(ji,jj,jk) * vmask(ji,jj,jk-1)227 zww(ji,jj,jk) = zww(ji,jj,jk) * wvmask(ji,jj,jk) 217 228 218 229 … … 301 312 zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj), 10._wp ) 302 313 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) & 303 & + zck * wslpiml(ji,jj) * zfk ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)314 & + zck * wslpiml(ji,jj) * zfk ) * wmask(ji,jj,jk) 304 315 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) & 305 & + zck * wslpjml(ji,jj) * zfk ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)316 & + zck * wslpjml(ji,jj) * zfk ) * wmask(ji,jj,jk) 306 317 307 318 !!gm modif to suppress omlmask.... (as in Griffies operator) … … 356 367 zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & 357 368 & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 358 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck * tmask(ji,jj,jk-1) * tmask(ji,jj,jk)359 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck * tmask(ji,jj,jk-1) * tmask(ji,jj,jk)369 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck * wmask(ji,jj,jk) 370 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck * wmask(ji,jj,jk) 360 371 END DO 361 372 END DO … … 423 434 vslp(ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 424 435 wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) & 425 & * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * 0.5436 & * wmask(ji,jj,jk) * 0.5 426 437 wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) & 427 & * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * 0.5438 & * wmask(ji,jj,jk) * 0.5 428 439 END DO 429 440 END DO … … 794 805 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj ) ) 795 806 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 796 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik)797 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik)807 wslpiml(ji,jj) = zai / ( zbi - zeps ) * wmask (ji,jj,ik) 808 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * wmask (ji,jj,ik) 798 809 END DO 799 810 END DO -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
r4990 r5098 98 98 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: emp_tot !: total E-P over ocean and ice [Kg/m2/s] 99 99 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fmmflx !: freshwater budget: freezing/melting [Kg/m2/s] 100 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rnf , rnf_b !: river runoff [Kg/m2/s] 100 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rnf , rnf_b !: river runoff [Kg/m2/s] 101 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf , fwfisf_b !: ice shelf melting [Kg/m2/s] 101 102 !! 102 103 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sbc_tsc, sbc_tsc_b !: sbc content trend [K.m/s] jpi,jpj,jpts … … 147 148 & sfx (jpi,jpj) , sfx_b(jpi,jpj) , emp_tot(jpi,jpj), fmmflx(jpi,jpj), STAT=ierr(2) ) 148 149 ! 149 ALLOCATE( rnf (jpi,jpj) , sbc_tsc (jpi,jpj,jpts) , qsr_hc (jpi,jpj,jpk) , &150 & rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) )150 ALLOCATE( fwfisf (jpi,jpj), rnf (jpi,jpj) , sbc_tsc (jpi,jpj,jpts) , qsr_hc (jpi,jpj,jpk) , & 151 & fwfisf_b(jpi,jpj), rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) ) 151 152 ! 152 153 ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) , & -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
r4990 r5098 8 8 !! 3.0 ! 2006-08 (G. Madec) Surface module 9 9 !! 3.2 ! 2009-07 (C. Talandier) emp mean s spread over erp area 10 !! 3.6 ! 2014-11 (P. Mathiot ) add ice shelf melting 10 11 !!---------------------------------------------------------------------- 11 12 … … 158 159 zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) ) 159 160 ! ! fwf global mean (excluding ocean to ice/snow exchanges) 160 z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - snwice_fmass(:,:) ) ) / area161 z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) - snwice_fmass(:,:) ) ) / area 161 162 ! 162 163 IF( z_fwf < 0._wp ) THEN ! spread out over >0 erp area to increase evaporation -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
r4990 r5098 7 7 !! History : 3.2 ! 2011-02 (C.Harris ) Original code isf cav 8 8 !! X.X ! 2006-02 (C. Wang ) Original code bg03 9 !! 3.4 ! 2013-03 (P. Mathiot) Merging 9 !! 3.4 ! 2013-03 (P. Mathiot) Merging + parametrization 10 10 !!---------------------------------------------------------------------- 11 11 … … 37 37 38 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_tsc_b, risf_tsc 39 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf_b, fwfisf !: evaporation damping [kg/m2/s] 40 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qisf !: net heat flux from ice shelf 39 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qisf !: net heat flux from ice shelf 41 40 REAL(wp), PUBLIC :: rn_hisf_tbl !: thickness of top boundary layer [m] 42 41 LOGICAL , PUBLIC :: ln_divisf !: flag to correct divergence … … 309 308 sbc_isf_alloc = 0 ! set to zero if no array to be allocated 310 309 IF( .NOT. ALLOCATED( qisf ) ) THEN 311 ALLOCATE( risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts) , & 312 & qisf(jpi,jpj) , fwfisf(jpi,jpj) , fwfisf_b(jpi,jpj) , & 313 & rhisf_tbl(jpi,jpj), r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj) , & 314 & ttbl(jpi,jpj) , stbl(jpi,jpj) , utbl(jpi,jpj) , & 315 & vtbl(jpi, jpj) , risfLeff(jpi,jpj) , rhisf_tbl_0(jpi,jpj), & 316 & ralpha(jpi,jpj) , misfkt(jpi,jpj) , misfkb(jpi,jpj) , & 310 ALLOCATE( risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj) , & 311 & rhisf_tbl(jpi,jpj) , r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj) , & 312 & ttbl(jpi,jpj) , stbl(jpi,jpj) , utbl(jpi,jpj) , & 313 & vtbl(jpi, jpj) , risfLeff(jpi,jpj) , rhisf_tbl_0(jpi,jpj), & 314 & ralpha(jpi,jpj) , misfkt(jpi,jpj) , misfkb(jpi,jpj) , & 317 315 & STAT= sbc_isf_alloc ) 318 316 ! -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
r4990 r5098 13 13 !! 3.4 ! 2011-11 (C. Harris) CICE added as an option 14 14 !! 3.5 ! 2012-11 (A. Coward, G. Madec) Rethink of heat, mass and salt surface fluxes 15 !! 3.6 ! 2014-11 (P. Mathiot, C. Harris) add ice shelves melting 15 16 !!---------------------------------------------------------------------- 16 17 … … 179 180 IF( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' ) 180 181 fwfisf (:,:) = 0.0_wp 182 fwfisf_b(:,:) = 0.0_wp 181 183 END IF 182 184 IF( nn_ice == 0 ) fr_i(:,:) = 0.e0 ! no ice in the domain, ice fraction is always zero -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
r4990 r5098 106 106 ENDIF 107 107 ! 108 zwi(:,:,:) = 0.e0 ; zwz(:,:,:) = 0.e0108 zwi(:,:,:) = 0.e0 ; 109 109 ! 110 110 ! ! =========== 111 111 DO jn = 1, kjpt ! tracer loop 112 112 ! ! =========== 113 ! 1. Bottom value : flux set to zero113 ! 1. Bottom and k=1 value : flux set to zero 114 114 ! ---------------------------------- 115 115 zwx(:,:,jpk) = 0.e0 ; zwz(:,:,jpk) = 0.e0 116 116 zwy(:,:,jpk) = 0.e0 ; zwi(:,:,jpk) = 0.e0 117 117 118 zwz(:,:,1 ) = 0._wp 118 119 ! 2. upstream advection with initial mass fluxes & intermediate update 119 120 ! -------------------------------------------------------------------- … … 134 135 135 136 ! upstream tracer flux in the k direction 137 ! Interior value 138 DO jk = 2, jpkm1 139 DO jj = 1, jpj 140 DO ji = 1, jpi 141 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 142 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 143 zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 144 END DO 145 END DO 146 END DO 136 147 ! Surface value 137 148 IF( lk_vvl ) THEN 138 149 DO jj = 1, jpj 139 150 DO ji = 1, jpi 140 zwz(ji,jj, mikt(ji,jj) ) = 0.e0 151 zwz(ji,jj, mikt(ji,jj) ) = 0.e0 ! volume variable 141 152 END DO 142 153 END DO … … 148 159 END DO 149 160 ENDIF 150 ! Interior value151 DO jj = 1, jpj152 DO ji = 1, jpi153 DO jk = mikt(ji,jj)+1, jpkm1154 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )155 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )156 zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) )157 END DO158 END DO159 END DO160 161 161 162 ! total advective trend … … 202 203 203 204 ! antidiffusive flux on k 204 zwz(:,:,1) = 0.e0 ! Surface value 205 ! 205 ! Interior value 206 DO jk = 2, jpkm1 207 DO jj = 1, jpj 208 DO ji = 1, jpi 209 zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk) 210 END DO 211 END DO 212 END DO 213 ! surface value 206 214 DO jj = 1, jpj 207 215 DO ji = 1, jpi 208 ik=mikt(ji,jj) 209 ! surface value 210 zwz(ji,jj,1:ik) = 0.e0 211 ! Interior value 212 DO jk = mikt(ji,jj)+1, jpkm1 213 zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk) 214 END DO 216 zwz(ji,jj,mikt(ji,jj)) = 0.e0 215 217 END DO 216 218 END DO … … 580 582 & paft * tmask + zbig * ( 1._wp - tmask ) ) 581 583 582 DO j j = 2, jpjm1583 DO j i = fs_2, fs_jpim1 ! vector opt.584 DO j k = mikt(ji,jj), jpkm1585 ikm1 = MAX(jk-1, mikt(ji,jj))584 DO jk = 1, jpkm1 585 DO jj = 2, jpjm1 586 DO ji = fs_2, fs_jpim1 ! vector opt. 587 ikm1 = MAX(jk-1,1) 586 588 z2dtt = p2dt(jk) 587 589 -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90
r4990 r5098 116 116 END DO 117 117 END DO 118 119 118 ! !== Laplacian ==! 120 119 ! … … 125 124 END DO 126 125 END DO 126 ! 127 127 IF( ln_zps ) THEN ! set gradient at partial step level (last ocean level) 128 128 DO jj = 1, jpjm1 … … 130 130 IF( mbku(ji,jj) == jk ) ztu(ji,jj,jk) = zeeu(ji,jj) * pgu(ji,jj,jn) 131 131 IF( mbkv(ji,jj) == jk ) ztv(ji,jj,jk) = zeev(ji,jj) * pgv(ji,jj,jn) 132 ! (ISH)133 IF( miku(ji,jj) == jk ) ztu(ji,jj,jk) = zeeu(ji,jj) * pgui(ji,jj,jn)134 IF( mikv(ji,jj) == jk ) ztu(ji,jj,jk) = zeev(ji,jj) * pgvi(ji,jj,jn)135 132 END DO 136 133 END DO 137 134 ENDIF 135 ! (ISH) 136 IF( ln_zps .AND. ln_isfcav ) THEN ! set gradient at partial step level (first ocean level in a cavity) 137 DO jj = 1, jpjm1 138 DO ji = 1, jpim1 139 IF( miku(ji,jj) == MAX(jk,2) ) ztu(ji,jj,jk) = zeeu(ji,jj) * pgui(ji,jj,jn) 140 IF( mikv(ji,jj) == MAX(jk,2) ) ztu(ji,jj,jk) = zeev(ji,jj) * pgvi(ji,jj,jn) 141 END DO 142 END DO 143 ENDIF 144 ! 138 145 DO jj = 2, jpjm1 ! Second derivative (divergence) time the eddy diffusivity coefficient 139 146 DO ji = fs_2, fs_jpim1 ! vector opt. -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
r4990 r5098 106 106 ! 107 107 INTEGER :: ji, jj, jk, jn ! dummy loop indices 108 INTEGER :: ikt 108 109 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars 109 110 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - … … 149 150 END DO 150 151 END DO 152 153 ! partial cell correction 151 154 IF( ln_zps ) THEN ! partial steps correction at the last ocean level 152 155 DO jj = 1, jpjm1 153 156 DO ji = 1, fs_jpim1 ! vector opt. 154 157 ! IF useless if zpshde defines pgu everywhere 155 IF (mbku(ji,jj) > 1) zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 156 IF (mbkv(ji,jj) > 1) zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 157 ! (ISF) 158 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 159 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 160 END DO 161 END DO 162 ENDIF 163 IF( ln_zps .AND. ln_isfcav ) THEN ! partial steps correction at the first wet level beneath a cavity 164 DO jj = 1, jpjm1 165 DO ji = 1, fs_jpim1 ! vector opt. 158 166 IF (miku(ji,jj) > 1) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 159 167 IF (mikv(ji,jj) > 1) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 160 168 END DO 161 169 END DO 162 END IF170 END IF 163 171 164 172 !!---------------------------------------------------------------------- 165 173 !! II - horizontal trend (full) 166 174 !!---------------------------------------------------------------------- 167 !CDIR PARALLEL DO PRIVATE( zdk1t ) 168 ! ! =============== 169 DO jj = 1, jpj ! Horizontal slab 170 ! ! =============== 171 DO ji = 1, jpi ! vector opt. 172 DO jk = mikt(ji,jj), jpkm1 173 ! 1. Vertical tracer gradient at level jk and jk+1 174 ! ------------------------------------------------ 175 ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 176 zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn) - ptb(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1) 177 ! 178 IF( jk == mikt(ji,jj) ) THEN ; zdkt(ji,jj,jk) = zdk1t(ji,jj,jk) 179 ELSE ; zdkt(ji,jj,jk) = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 180 ENDIF 175 !!!!!!!!!!CDIR PARALLEL DO PRIVATE( zdk1t ) 176 ! 1. Vertical tracer gradient at level jk and jk+1 177 ! ------------------------------------------------ 178 ! 179 ! interior value 180 DO jk = 2, jpkm1 181 DO jj = 1, jpj 182 DO ji = 1, jpi ! vector opt. 183 zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn ) - ptb(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1) 184 ! 185 zdkt(ji,jj,jk) = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn ) ) * wmask(ji,jj,jk) 181 186 END DO 182 187 END DO 183 188 END DO 189 ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 190 zdk1t(:,:,1) = ( ptb(:,:,1,jn ) - ptb(:,:,2,jn) ) * wmask(:,:,2) 191 zdkt (:,:,1) = zdk1t(:,:,1) 192 IF ( ln_isfcav ) THEN 193 DO jj = 1, jpj 194 DO ji = 1, jpi ! vector opt. 195 ikt = mikt(ji,jj) ! surface level 196 zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1) 197 zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt) 198 END DO 199 END DO 200 END IF 184 201 185 202 ! 2. Horizontal fluxes 186 203 ! -------------------- 187 DO j j = 1 , jpjm1188 DO j i = 1, fs_jpim1 ! vector opt.189 DO j k = mikt(ji,jj), jpkm1204 DO jk = 1, jpkm1 205 DO jj = 1 , jpjm1 206 DO ji = 1, fs_jpim1 ! vector opt. 190 207 zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) 191 208 zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) … … 208 225 END DO 209 226 END DO 210 END DO211 227 212 228 ! II.4 Second derivative (divergence) and add to the general trend 213 229 ! ---------------------------------------------------------------- 214 DO jj = 2 , jpjm1 215 DO ji = fs_2, fs_jpim1 ! vector opt. 216 DO jk = mikt(ji,jj), jpkm1 217 zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 230 DO jj = 2 , jpjm1 231 DO ji = fs_2, fs_jpim1 ! vector opt. 232 zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 218 233 ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) 219 234 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra … … 278 293 DO jj = 2, jpjm1 279 294 DO ji = fs_2, fs_jpim1 ! vector opt. 280 zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)295 zcoef0 = - fsahtw(ji,jj,jk) * wmask(ji,jj,jk) 281 296 ! 282 297 zmsku = 1./MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90
r4990 r5098 102 102 END DO 103 103 END DO 104 IF( ln_zps ) THEN ! set gradient at partial step level 104 IF( ln_zps ) THEN ! set gradient at partial step level for the last ocean cell 105 105 DO jj = 1, jpjm1 106 106 DO ji = 1, fs_jpim1 ! vector opt. … … 116 116 ztv(ji,jj,jk) = zabe2 * pgv(ji,jj,jn) 117 117 ENDIF 118 119 ! (ISH) 118 END DO 119 END DO 120 ENDIF 121 ! (ISH) 122 IF( ln_zps .AND. ln_isfcav ) THEN ! set gradient at partial step level for the first ocean cell 123 ! into a cavity 124 DO jj = 1, jpjm1 125 DO ji = 1, fs_jpim1 ! vector opt. 120 126 ! ice shelf level level MAX(2,jk) => only where ice shelf 121 127 iku = miku(ji,jj) -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
r4990 r5098 9 9 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps 10 10 !! - ! 2010-09 (C. Ethe, G. Madec) Merge TRA-TRC 11 !! 3.6 ! 2014-11 (P. Mathiot) isf melting forcing 11 12 !!---------------------------------------------------------------------- 12 13 -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90
r4990 r5098 122 122 DO jj=1, jpj 123 123 DO ji=1, jpi 124 zwt(ji,jj,1 :mikt(ji,jj)) = 0._wp124 zwt(ji,jj,1) = 0._wp 125 125 END DO 126 126 END DO … … 184 184 DO jj = 2, jpjm1 185 185 DO ji = fs_2, fs_jpim1 186 zwt(ji,jj,1:mikt(ji,jj)) = zwd(ji,jj,1:mikt(ji,jj)) 187 DO jk = mikt(ji,jj)+1, jpkm1 186 zwt(ji,jj,1) = zwd(ji,jj,1) 187 END DO 188 END DO 189 DO jk = 2, jpkm1 190 DO jj = 2, jpjm1 191 DO ji = fs_2, fs_jpim1 188 192 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 189 193 END DO … … 196 200 DO jj = 2, jpjm1 197 201 DO ji = fs_2, fs_jpim1 198 ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,mikt(ji,jj)) 199 ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,mikt(ji,jj)) 200 pta(ji,jj,mikt(ji,jj),jn) = ze3tb * ptb(ji,jj,mikt(ji,jj),jn) & 201 & + p2dt(mikt(ji,jj)) * ze3tn * pta(ji,jj,mikt(ji,jj),jn) 202 DO jk = mikt(ji,jj)+1, jpkm1 202 ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1) 203 ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1) 204 pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) & 205 & + p2dt(1) * ze3tn * pta(ji,jj,1,jn) 206 END DO 207 END DO 208 DO jk = 2, jpkm1 209 DO jj = 2, jpjm1 210 DO ji = fs_2, fs_jpim1 203 211 ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,jk) 204 212 ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t (ji,jj,jk) … … 213 221 DO ji = fs_2, fs_jpim1 214 222 pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 215 DO jk = jpk-2, mikt(ji,jj), -1 223 END DO 224 END DO 225 DO jk = jpk-2, 1, -1 226 DO jj = 2, jpjm1 227 DO ji = fs_2, fs_jpim1 216 228 pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) & 217 229 & / zwt(ji,jj,jk) * tmask(ji,jj,jk) -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
r4990 r5098 8 8 !! - ! 2004-03 (C. Ethe) adapted for passive tracers 9 9 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA 10 !! 3.6 ! 2014-11 (P. Mathiot) Add zps_hde_isf (needed to open a cavity) 10 11 !!====================================================================== 11 12 … … 27 28 PRIVATE 28 29 29 PUBLIC zps_hde ! routine called by step.F90 30 PUBLIC zps_hde ! routine called by step.F90 31 PUBLIC zps_hde_isf ! routine called by step.F90 30 32 31 33 !! * Substitutions … … 40 42 41 43 SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv, & 42 & prd, pgru, pgrv, pmru, pmrv, pgzu, pgzv, pge3ru, pge3rv, & 43 & sgtu, sgtv, sgru, sgrv, smru, smrv, sgzu, sgzv, sge3ru, sge3rv ) 44 & prd, pgru, pgrv ) 44 45 !!---------------------------------------------------------------------- 45 46 !! *** ROUTINE zps_hde *** … … 82 83 !! 83 84 !! ** Action : compute for top and bottom interfaces 84 !! - pgtu, pgtv, sgtu, sgtv: horizontal gradient of tracer at u- & v-points 85 !! - pgru, pgrv, sgru, sgtv: horizontal gradient of rho (if present) at u- & v-points 86 !! - pmru, pmrv, smru, smrv: horizontal sum of rho at u- & v- point (used in dynhpg with vvl) 87 !! - pgzu, pgzv, sgzu, sgzv: horizontal gradient of z at u- and v- point (used in dynhpg with vvl) 88 !! - pge3ru, pge3rv, sge3ru, sge3rv: horizontal gradient of rho weighted by local e3w at u- & v-points 85 !! - pgtu, pgtv: horizontal gradient of tracer at u- & v-points 86 !! - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points 89 87 !!---------------------------------------------------------------------- 90 88 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 92 90 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pta ! 4D tracers fields 93 91 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 94 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: sgtu, sgtv ! hor. grad. of stra at u- & v-pts (ISF)95 92 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 96 93 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 97 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pmru, pmrv ! hor. sum of prd at u- & v-pts (bottom)98 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgzu, pgzv ! hor. grad of z at u- & v-pts (bottom)99 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pge3ru, pge3rv ! hor. grad of prd weighted by local e3w at u- & v-pts (bottom)100 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: sgru, sgrv ! hor. grad of prd at u- & v-pts (top)101 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: smru, smrv ! hor. sum of prd at u- & v-pts (top)102 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: sgzu, sgzv ! hor. grad of z at u- & v-pts (top)103 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: sge3ru, sge3rv ! hor. grad of prd weighted by local e3w at u- & v-pts (top)104 94 ! 105 95 INTEGER :: ji, jj, jn ! Dummy loop indices … … 113 103 ! 114 104 pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 115 sgtu(:,:,:)=0.0_wp ; sgtv(:,:,:)=0.0_wp ; 105 zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 106 zhi (:,: )=0.0_wp ; zhj (:,: )=0.0_wp ; 107 ! 108 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! 109 ! 110 DO jj = 1, jpjm1 111 DO ji = 1, jpim1 112 iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 113 ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 ) ! if level first is a p-step, ik.m1=1 114 ! (ISF) case partial step top and bottom in adjacent cell in vertical 115 ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 116 ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 117 ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 118 ze3wu = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 119 ze3wv = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 120 ! 121 ! i- direction 122 IF( ze3wu >= 0._wp ) THEN ! case 1 123 zmaxu = ze3wu / fse3w(ji+1,jj,iku) 124 ! interpolated values of tracers 125 zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 126 ! gradient of tracers 127 pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 128 ELSE ! case 2 129 zmaxu = -ze3wu / fse3w(ji,jj,iku) 130 ! interpolated values of tracers 131 zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 132 ! gradient of tracers 133 pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 134 ENDIF 135 ! 136 ! j- direction 137 IF( ze3wv >= 0._wp ) THEN ! case 1 138 zmaxv = ze3wv / fse3w(ji,jj+1,ikv) 139 ! interpolated values of tracers 140 ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 141 ! gradient of tracers 142 pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 143 ELSE ! case 2 144 zmaxv = -ze3wv / fse3w(ji,jj,ikv) 145 ! interpolated values of tracers 146 ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 147 ! gradient of tracers 148 pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 149 ENDIF 150 END DO 151 END DO 152 CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. ) ! Lateral boundary cond. 153 ! 154 END DO 155 156 ! horizontal derivative of density anomalies (rd) 157 IF( PRESENT( prd ) ) THEN ! depth of the partial step level 158 pgru(:,:)=0.0_wp ; pgrv(:,:)=0.0_wp ; 159 DO jj = 1, jpjm1 160 DO ji = 1, jpim1 161 iku = mbku(ji,jj) 162 ikv = mbkv(ji,jj) 163 ze3wu = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 164 ze3wv = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 165 166 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu ! i-direction: case 1 167 ELSE ; zhi(ji,jj) = fsdept(ji ,jj,iku) + ze3wu ! - - case 2 168 ENDIF 169 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = fsdept(ji,jj+1,ikv) - ze3wv ! j-direction: case 1 170 ELSE ; zhj(ji,jj) = fsdept(ji,jj ,ikv) + ze3wv ! - - case 2 171 ENDIF 172 END DO 173 END DO 174 175 ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 176 ! step and store it in zri, zrj for each case 177 CALL eos( zti, zhi, zri ) 178 CALL eos( ztj, zhj, zrj ) 179 180 ! Gradient of density at the last level 181 DO jj = 1, jpjm1 182 DO ji = 1, jpim1 183 iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 184 ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 ) ! last and before last ocean level at u- & v-points 185 ze3wu = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 186 ze3wv = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 187 IF( ze3wu >= 0._wp ) THEN ; pgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji ,jj) - prd(ji,jj,iku) ) ! i: 1 188 ELSE ; pgru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) ) ! i: 2 189 ENDIF 190 IF( ze3wv >= 0._wp ) THEN ; pgrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 191 ELSE ; pgrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) ) ! j: 2 192 ENDIF 193 END DO 194 END DO 195 CALL lbc_lnk( pgru , 'U', -1. ) ; CALL lbc_lnk( pgrv , 'V', -1. ) ! Lateral boundary conditions 196 ! 197 END IF 198 ! 199 IF( nn_timing == 1 ) CALL timing_stop( 'zps_hde') 200 ! 201 END SUBROUTINE zps_hde 202 ! 203 SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, & 204 & prd, pgru, pgrv, pmru, pmrv, pgzu, pgzv, pge3ru, pge3rv, & 205 & pgtui, pgtvi, pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi ) 206 !!---------------------------------------------------------------------- 207 !! *** ROUTINE zps_hde *** 208 !! 209 !! ** Purpose : Compute the horizontal derivative of T, S and rho 210 !! at u- and v-points with a linear interpolation for z-coordinate 211 !! with partial steps. 212 !! 213 !! ** Method : In z-coord with partial steps, scale factors on last 214 !! levels are different for each grid point, so that T, S and rd 215 !! points are not at the same depth as in z-coord. To have horizontal 216 !! gradients again, we interpolate T and S at the good depth : 217 !! Linear interpolation of T, S 218 !! Computation of di(tb) and dj(tb) by vertical interpolation: 219 !! di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~ 220 !! dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~ 221 !! This formulation computes the two cases: 222 !! CASE 1 CASE 2 223 !! k-1 ___ ___________ k-1 ___ ___________ 224 !! Ti T~ T~ Ti+1 225 !! _____ _____ 226 !! k | |Ti+1 k Ti | | 227 !! | |____ ____| | 228 !! ___ | | | ___ | | | 229 !! 230 !! case 1-> e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then 231 !! t~ = t(i+1,j ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1) 232 !! ( t~ = t(i ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1) ) 233 !! or 234 !! case 2-> e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then 235 !! t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i ) 236 !! ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) ) 237 !! Idem for di(s) and dj(s) 238 !! 239 !! For rho, we call eos which will compute rd~(t~,s~) at the right 240 !! depth zh from interpolated T and S for the different formulations 241 !! of the equation of state (eos). 242 !! Gradient formulation for rho : 243 !! di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~ 244 !! 245 !! ** Action : compute for top and bottom interfaces 246 !! - pgtu, pgtv, pgtui, pgtvi: horizontal gradient of tracer at u- & v-points 247 !! - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points 248 !! - pmru, pmrv, pmrui, pmrvi: horizontal sum of rho at u- & v- point (used in dynhpg with vvl) 249 !! - pgzu, pgzv, pgzui, pgzvi: horizontal gradient of z at u- and v- point (used in dynhpg with vvl) 250 !! - pge3ru, pge3rv, pge3rui, pge3rvi: horizontal gradient of rho weighted by local e3w at u- & v-points 251 !!---------------------------------------------------------------------- 252 INTEGER , INTENT(in ) :: kt ! ocean time-step index 253 INTEGER , INTENT(in ) :: kjpt ! number of tracers 254 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pta ! 4D tracers fields 255 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 256 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF) 257 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 258 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 259 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pmru, pmrv ! hor. sum of prd at u- & v-pts (bottom) 260 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgzu, pgzv ! hor. grad of z at u- & v-pts (bottom) 261 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pge3ru, pge3rv ! hor. grad of prd weighted by local e3w at u- & v-pts (bottom) 262 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top) 263 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pmrui, pmrvi ! hor. sum of prd at u- & v-pts (top) 264 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgzui, pgzvi ! hor. grad of z at u- & v-pts (top) 265 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pge3rui, pge3rvi ! hor. grad of prd weighted by local e3w at u- & v-pts (top) 266 ! 267 INTEGER :: ji, jj, jn ! Dummy loop indices 268 INTEGER :: iku, ikv, ikum1, ikvm1,ikup1, ikvp1 ! partial step level (ocean bottom level) at u- and v-points 269 REAL(wp) :: ze3wu, ze3wv, zmaxu, zmaxv, zdzwu, zdzwv, zdzwuip1, zdzwvjp1 ! temporary scalars 270 REAL(wp), DIMENSION(jpi,jpj) :: zri, zrj, zhi, zhj ! NB: 3rd dim=1 to use eos 271 REAL(wp), DIMENSION(jpi,jpj,kjpt) :: zti, ztj ! 272 !!---------------------------------------------------------------------- 273 ! 274 IF( nn_timing == 1 ) CALL timing_start( 'zps_hde_isf') 275 ! 276 pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 277 pgtui(:,:,:)=0.0_wp ; pgtvi(:,:,:)=0.0_wp ; 116 278 zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 117 279 zhi (:,: )=0.0_wp ; zhj (:,: )=0.0_wp ; … … 256 418 zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) ) 257 419 ! gradient of tracers 258 sgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )420 pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 259 421 ELSE ! case 2 260 422 zmaxu = - ze3wu / fse3w(ji,jj,iku+1) … … 262 424 zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) ) 263 425 ! gradient of tracers 264 sgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )426 pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 265 427 ENDIF 266 428 ! … … 271 433 ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) ) 272 434 ! gradient of tracers 273 sgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )435 pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 274 436 ELSE ! case 2 275 437 zmaxv = - ze3wv / fse3w(ji,jj,ikv+1) … … 277 439 ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) ) 278 440 ! gradient of tracers 279 sgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )441 pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 280 442 ENDIF 281 443 END DO!! 282 444 END DO!! 283 CALL lbc_lnk( sgtu(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( sgtv(:,:,jn), 'V', -1. ) ! Lateral boundary cond.445 CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. ) ! Lateral boundary cond. 284 446 ! 285 447 END DO … … 287 449 ! horizontal derivative of density anomalies (rd) 288 450 IF( PRESENT( prd ) ) THEN ! depth of the partial step level 289 sgru(:,:) =0.0_wp ; sgrv(:,:) =0.0_wp ;290 sgzu(:,:) =0.0_wp ; sgzv(:,:) =0.0_wp ;291 smru(:,:) =0.0_wp ; smru(:,:) =0.0_wp ;292 sge3ru(:,:)=0.0_wp ; sge3rv(:,:)=0.0_wp ;451 pgrui(:,:) =0.0_wp ; pgrvi(:,:) =0.0_wp ; 452 pgzui(:,:) =0.0_wp ; pgzvi(:,:) =0.0_wp ; 453 pmrui(:,:) =0.0_wp ; pmrui(:,:) =0.0_wp ; 454 pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ; 293 455 294 456 DO jj = 1, jpjm1 … … 321 483 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)) 322 484 IF( ze3wu >= 0._wp ) THEN 323 sgzu(ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku)324 sgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji,jj) - prd(ji,jj,iku) ) ! i: 1325 smru(ji,jj) = umask(ji,jj,iku) * ( zri(ji,jj) + prd(ji,jj,iku) ) ! i: 1326 sge3ru(ji,jj) = umask(ji,jj,iku+1) &485 pgzui (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku) 486 pgrui (ji,jj) = umask(ji,jj,iku) * ( zri(ji,jj) - prd(ji,jj,iku) ) ! i: 1 487 pmrui (ji,jj) = umask(ji,jj,iku) * ( zri(ji,jj) + prd(ji,jj,iku) ) ! i: 1 488 pge3rui(ji,jj) = umask(ji,jj,iku+1) & 327 489 * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj ) + prd(ji+1,jj,iku+1) + 2._wp) & 328 490 - fse3w(ji ,jj,iku+1) * (prd(ji,jj,iku) + prd(ji ,jj,iku+1) + 2._wp) ) ! i: 1 329 491 ELSE 330 sgzu(ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu)331 sgru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) ) ! i: 2332 smru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) ) ! i: 2333 sge3ru(ji,jj) = umask(ji,jj,iku+1) &492 pgzui (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu) 493 pgrui (ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) ) ! i: 2 494 pmrui (ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) ) ! i: 2 495 pge3rui(ji,jj) = umask(ji,jj,iku+1) & 334 496 * ( fse3w(ji+1,jj,iku+1) * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp) & 335 497 -(fse3w(ji ,jj,iku+1) + ze3wu) * (zri(ji,jj ) + prd(ji ,jj,iku+1) + 2._wp) ) ! i: 2 336 498 ENDIF 337 499 IF( ze3wv >= 0._wp ) THEN 338 sgzv(ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv)339 sgrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1340 smrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) + prd(ji,jj,ikv) ) ! j: 1341 sge3rv(ji,jj) = vmask(ji,jj,ikv+1) &500 pgzvi (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv) 501 pgrvi (ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 502 pmrvi (ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) + prd(ji,jj,ikv) ) ! j: 1 503 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1) & 342 504 * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj ) + prd(ji,jj+1,ikv+1) + 2._wp) & 343 505 - fse3w(ji,jj ,ikv+1) * ( prd(ji,jj,ikv) + prd(ji,jj ,ikv+1) + 2._wp) ) ! j: 1 344 506 ! + 2 due to the formulation in density and not in anomalie in hpg sco 345 507 ELSE 346 sgzv(ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv)347 sgrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) ) ! j: 2348 smrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) ) ! j: 2349 sge3rv(ji,jj) = vmask(ji,jj,ikv+1) &508 pgzvi (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv) 509 pgrvi (ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) ) ! j: 2 510 pmrvi (ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) ) ! j: 2 511 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1) & 350 512 * ( fse3w(ji,jj+1,ikv+1) * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) & 351 513 -(fse3w(ji,jj ,ikv+1) + ze3wv) * ( zrj(ji,jj ) + prd(ji,jj ,ikv+1) + 2._wp) ) ! j: 2 … … 353 515 END DO 354 516 END DO 355 CALL lbc_lnk( sgru , 'U', -1. ) ; CALL lbc_lnk( sgrv, 'V', -1. ) ! Lateral boundary conditions356 CALL lbc_lnk( smru , 'U', 1. ) ; CALL lbc_lnk( smrv, 'V', 1. ) ! Lateral boundary conditions357 CALL lbc_lnk( sgzu , 'U', -1. ) ; CALL lbc_lnk( sgzv, 'V', -1. ) ! Lateral boundary conditions358 CALL lbc_lnk( sge3ru , 'U', -1. ) ; CALL lbc_lnk( sge3rv, 'V', -1. ) ! Lateral boundary conditions517 CALL lbc_lnk( pgrui , 'U', -1. ) ; CALL lbc_lnk( pgrvi , 'V', -1. ) ! Lateral boundary conditions 518 CALL lbc_lnk( pmrui , 'U', 1. ) ; CALL lbc_lnk( pmrvi , 'V', 1. ) ! Lateral boundary conditions 519 CALL lbc_lnk( pgzui , 'U', -1. ) ; CALL lbc_lnk( pgzvi , 'V', -1. ) ! Lateral boundary conditions 520 CALL lbc_lnk( pge3rui , 'U', -1. ) ; CALL lbc_lnk( pge3rvi , 'V', -1. ) ! Lateral boundary conditions 359 521 ! 360 522 END IF 361 523 ! 362 IF( nn_timing == 1 ) CALL timing_stop( 'zps_hde') 363 ! 364 END SUBROUTINE zps_hde 365 524 IF( nn_timing == 1 ) CALL timing_stop( 'zps_hde_isf') 525 ! 526 END SUBROUTINE zps_hde_isf 366 527 !!====================================================================== 367 528 END MODULE zpshde -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90
r4990 r5098 156 156 END DO 157 157 ! mask zmsk in order to have avt and avs masked 158 zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk)158 zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk) 159 159 160 160 … … 191 191 avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), & 192 192 & avt(ji,jj,jk), avt(ji+1,jj,jk), & 193 & avs(ji,jj,jk), avs(ji+1,jj,jk) ) * umask(ji,jj,jk)193 & avs(ji,jj,jk), avs(ji+1,jj,jk) ) * wumask(ji,jj,jk) 194 194 avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), & 195 195 & avt(ji,jj,jk), avt(ji,jj+1,jk), & 196 & avs(ji,jj,jk), avs(ji,jj+1,jk) ) * vmask(ji,jj,jk)196 & avs(ji,jj,jk), avs(ji,jj+1,jk) ) * wvmask(ji,jj,jk) 197 197 END DO 198 198 END DO … … 255 255 IF( zdf_ddm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_ddm_init : unable to allocate arrays' ) 256 256 ! ! initialization to masked Kz 257 avs(:,:,:) = rn_avt0 * tmask(:,:,:)257 avs(:,:,:) = rn_avt0 * wmask(:,:,:) 258 258 ! 259 259 END SUBROUTINE zdf_ddm_init -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r4990 r5098 26 26 !! ! + cleaning of the parameters + bugs correction 27 27 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 28 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 28 29 !!---------------------------------------------------------------------- 29 30 #if defined key_zdftke || defined key_esopa … … 236 237 zfact3 = 0.5_wp * rn_ediss 237 238 ! 239 ! 238 240 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 239 241 ! ! Surface boundary condition on tke 240 242 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 243 DO jj = 2, jpjm1 ! en(mikt(ji,jj)) = rn_emin 244 DO ji = fs_2, fs_jpim1 ! vector opt. 245 en(ji,jj,mikt(ji,jj))=rn_emin 246 END DO 247 END DO 241 248 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 242 249 DO ji = fs_2, fs_jpim1 ! vector opt. 243 IF (mikt(ji,jj) .GT. 1) THEN 244 en(ji,jj,mikt(ji,jj))=rn_emin 245 ELSE 246 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 247 END IF 250 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 248 251 END DO 249 252 END DO … … 301 304 END DO 302 305 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 306 !CDIR NOVERRCHK 303 307 DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en 304 DO jj = 2, jpjm1 308 !CDIR NOVERRCHK 309 DO jj = 2, jpjm1 310 !CDIR NOVERRCHK 305 311 DO ji = fs_2, fs_jpim1 ! vector opt. 306 312 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift … … 309 315 zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 310 316 ! ! TKE Langmuir circulation source term 311 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk)317 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 312 318 END DO 313 319 END DO … … 338 344 END DO 339 345 ! 340 DO j j = 2, jpjm1341 DO j i = fs_2, fs_jpim1 ! vector opt.342 DO j k = mikt(ji,jj)+1, jpkm1 !* Matrix and right hand side in en346 DO jk = 2, jpkm1 !* Matrix and right hand side in en 347 DO jj = 2, jpjm1 348 DO ji = fs_2, fs_jpim1 ! vector opt. 343 349 zcof = zfact1 * tmask(ji,jj,jk) 344 350 zzd_up = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & ! upper diagonal … … 357 363 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) & 358 364 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) & 359 & * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 360 END DO 361 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 362 DO jk = mikt(ji,jj)+2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 365 & * wmask(ji,jj,jk) 366 END DO 367 END DO 368 END DO 369 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 370 DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 371 DO jj = 2, jpjm1 372 DO ji = fs_2, fs_jpim1 ! vector opt. 363 373 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 364 374 END DO 365 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 366 zd_lw(ji,jj,mikt(ji,jj)+1) = en(ji,jj,mikt(ji,jj)+1) - zd_lw(ji,jj,mikt(ji,jj)+1) * en(ji,jj,mikt(ji,jj)) ! Surface boudary conditions on tke 367 ! 368 DO jk = mikt(ji,jj)+2, jpkm1 375 END DO 376 END DO 377 ! 378 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 379 DO jj = 2, jpjm1 380 DO ji = fs_2, fs_jpim1 ! vector opt. 381 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 382 END DO 383 END DO 384 DO jk = 3, jpkm1 385 DO jj = 2, jpjm1 386 DO ji = fs_2, fs_jpim1 ! vector opt. 369 387 zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1) 370 388 END DO 371 ! 372 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 389 END DO 390 END DO 391 ! 392 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 393 DO jj = 2, jpjm1 394 DO ji = fs_2, fs_jpim1 ! vector opt. 373 395 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 374 ! 375 DO jk = jpk-2, mikt(ji,jj)+1, -1 396 END DO 397 END DO 398 DO jk = jpk-2, 2, -1 399 DO jj = 2, jpjm1 400 DO ji = fs_2, fs_jpim1 ! vector opt. 376 401 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 377 402 END DO 378 ! 379 DO jk = mikt(ji,jj), jpkm1 ! set the minimum value of tke 380 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk) 403 END DO 404 END DO 405 DO jk = 2, jpkm1 ! set the minimum value of tke 406 DO jj = 2, jpjm1 407 DO ji = fs_2, fs_jpim1 ! vector opt. 408 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 381 409 END DO 382 410 END DO … … 391 419 DO ji = fs_2, fs_jpim1 ! vector opt. 392 420 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 393 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,1)421 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 394 422 END DO 395 423 END DO … … 400 428 jk = nmln(ji,jj) 401 429 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 402 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,1)430 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 403 431 END DO 404 432 END DO … … 416 444 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 417 445 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 418 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * tmask(ji,jj,1)446 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 419 447 END DO 420 448 END DO … … 484 512 ! !* Buoyancy length scale: l=sqrt(2*e/n**2) 485 513 ! 514 ! initialisation of interior minimum value (avoid a 2d loop with mikt) 515 zmxlm(:,:,:) = rmxl_min 516 zmxld(:,:,:) = rmxl_min 517 ! 486 518 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 487 519 DO jj = 2, jpjm1 488 520 DO ji = fs_2, fs_jpim1 489 IF (mikt(ji,jj) .GT. 1) THEN 490 zmxlm(ji,jj,mikt(ji,jj)) = rmxl_min 491 ELSE 492 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 493 zmxlm(ji,jj,mikt(ji,jj)) = MAX( rn_mxl0, zraug * taum(ji,jj) ) 494 END IF 521 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 522 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 495 523 END DO 496 524 END DO 497 525 ELSE 498 DO jj = 2, jpjm1 499 DO ji = fs_2, fs_jpim1 ! surface set to the minimum value 500 zmxlm(ji,jj,mikt(ji,jj)) = MAX( tmask(ji,jj,1) * rn_mxl0, rmxl_min) 501 END DO 502 END DO 526 zmxlm(:,:,1) = rn_mxl0 503 527 ENDIF 504 zmxlm(:,:,jpk) = rmxl_min ! last level set to the interior minium value 505 ! 506 !CDIR NOVERRCHK 507 DO jj = 2, jpjm1 508 !CDIR NOVERRCHK 509 DO ji = fs_2, fs_jpim1 ! vector opt. 510 !CDIR NOVERRCHK 511 DO jk = mikt(ji,jj)+1, jpkm1 ! interior value : l=sqrt(2*e/n^2) 528 ! 529 !CDIR NOVERRCHK 530 DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2) 531 !CDIR NOVERRCHK 532 DO jj = 2, jpjm1 533 !CDIR NOVERRCHK 534 DO ji = fs_2, fs_jpim1 ! vector opt. 512 535 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 513 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 514 END DO 515 zmxld(ji,jj,mikt(ji,jj)) = zmxlm(ji,jj,mikt(ji,jj)) ! surface set to the minimum value 536 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 537 END DO 516 538 END DO 517 539 END DO … … 519 541 ! !* Physical limits for the mixing length 520 542 ! 521 zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the zmxlm value543 zmxld(:,:,1 ) = zmxlm(:,:,1) ! surface set to the minimum value 522 544 zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value 523 545 ! 524 546 SELECT CASE ( nn_mxl ) 525 547 ! 548 ! where wmask = 0 set zmxlm == fse3w 526 549 CASE ( 0 ) ! bounded by the distance to surface and bottom 527 DO j j = 2, jpjm1528 DO j i = fs_2, fs_jpim1 ! vector opt.529 DO j k = mikt(ji,jj)+1, jpkm1550 DO jk = 2, jpkm1 551 DO jj = 2, jpjm1 552 DO ji = fs_2, fs_jpim1 ! vector opt. 530 553 zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 531 554 & fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 532 zmxlm(ji,jj,jk) = zemxl 533 zmxld(ji,jj,jk) = zemxl 555 ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 556 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 557 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 534 558 END DO 535 559 END DO … … 537 561 ! 538 562 CASE ( 1 ) ! bounded by the vertical scale factor 539 DO j j = 2, jpjm1540 DO j i = fs_2, fs_jpim1 ! vector opt.541 DO j k = mikt(ji,jj)+1, jpkm1563 DO jk = 2, jpkm1 564 DO jj = 2, jpjm1 565 DO ji = fs_2, fs_jpim1 ! vector opt. 542 566 zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 543 567 zmxlm(ji,jj,jk) = zemxl … … 548 572 ! 549 573 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 550 DO j j = 2, jpjm1551 DO j i = fs_2, fs_jpim1 ! vector opt.552 DO j k = mikt(ji,jj)+1, jpkm1 ! from the surface to the bottom :574 DO jk = 2, jpkm1 ! from the surface to the bottom : 575 DO jj = 2, jpjm1 576 DO ji = fs_2, fs_jpim1 ! vector opt. 553 577 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 554 578 END DO 555 DO jk = jpkm1, mikt(ji,jj)+1, -1 ! from the bottom to the surface : 579 END DO 580 END DO 581 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : 582 DO jj = 2, jpjm1 583 DO ji = fs_2, fs_jpim1 ! vector opt. 556 584 zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 557 585 zmxlm(ji,jj,jk) = zemxl … … 562 590 ! 563 591 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 564 DO j j = 2, jpjm1565 DO j i = fs_2, fs_jpim1 ! vector opt.566 DO j k = mikt(ji,jj)+1, jpkm1 ! from the surface to the bottom : lup592 DO jk = 2, jpkm1 ! from the surface to the bottom : lup 593 DO jj = 2, jpjm1 594 DO ji = fs_2, fs_jpim1 ! vector opt. 567 595 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 568 596 END DO 569 DO jk = jpkm1, mikt(ji,jj)+1, -1 ! from the bottom to the surface : ldown 597 END DO 598 END DO 599 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown 600 DO jj = 2, jpjm1 601 DO ji = fs_2, fs_jpim1 ! vector opt. 570 602 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 571 603 END DO … … 604 636 zsqen = SQRT( en(ji,jj,jk) ) 605 637 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen 606 avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)607 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)638 avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 639 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 608 640 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 609 641 END DO … … 612 644 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 613 645 ! 614 DO jj = 2, jpjm1 615 DO ji = fs_2, fs_jpim1 ! vector opt. 616 DO jk = miku(ji,jj)+1, jpkm1 !* vertical eddy viscosity at u- and v-points 617 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk) 618 END DO 619 DO jk = mikv(ji,jj)+1, jpkm1 620 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk) 646 DO jk = 2, jpkm1 !* vertical eddy viscosity at wu- and wv-points 647 DO jj = 2, jpjm1 648 DO ji = fs_2, fs_jpim1 ! vector opt. 649 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * wumask(ji,jj,jk) 650 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * wvmask(ji,jj,jk) 621 651 END DO 622 652 END DO … … 625 655 ! 626 656 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 627 DO j j = 2, jpjm1628 DO j i = fs_2, fs_jpim1 ! vector opt.629 DO j k = mikt(ji,jj)+1, jpkm1657 DO jk = 2, jpkm1 658 DO jj = 2, jpjm1 659 DO ji = fs_2, fs_jpim1 ! vector opt. 630 660 zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 631 661 ! ! shear … … 639 669 !!gm and even better with the use of the "true" ri_crit=0.22222... (this change the results!) 640 670 !!gm zpdlr = MAX( 0.1_wp, ri_crit / MAX( ri_crit , zri ) ) 641 avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)671 avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 642 672 # if defined key_c1d 643 e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number644 e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk) ! c1d config. : save Ri673 e_pdl(ji,jj,jk) = zpdlr * wmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number 674 e_ric(ji,jj,jk) = zri * wmask(ji,jj,jk) ! c1d config. : save Ri 645 675 # endif 646 676 END DO … … 749 779 ! !* set vertical eddy coef. to the background value 750 780 DO jk = 1, jpk 751 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)752 avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)753 avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)754 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)781 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 782 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 783 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 784 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 755 785 END DO 756 786 dissl(:,:,:) = 1.e-12_wp … … 808 838 en(:,:,:) = rn_emin * tmask(:,:,:) 809 839 DO jk = 1, jpk ! set the Kz to the background value 810 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)811 avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)812 avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)813 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)840 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 841 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 842 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 843 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 814 844 END DO 815 845 ENDIF -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90
r5021 r5098 126 126 zkz(:,:) = 0.e0 !* Associated potential energy consummed over the whole water column 127 127 DO jk = 2, jpkm1 128 zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)128 zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk) * wmask(:,:,jk) 129 129 END DO 130 130 … … 135 135 END DO 136 136 137 DO j j = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx138 DO j i = 1, jpi139 DO j k = mikt(ji,jj)+1, jpkm1 !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s140 zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) !kz max = 300 cm2/s137 DO jk = 2, jpkm1 !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s 138 DO jj = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 139 DO ji = 1, jpi 140 zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk) !kz max = 300 cm2/s 141 141 END DO 142 142 END DO … … 166 166 ! ! Update mixing coefs ! 167 167 ! ! ----------------------- ! 168 DO j j = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx169 DO j i = 1, jpi170 DO j k = mikt(ji,jj)+1, jpkm1 !* update momentum & tracer diffusivity with tidal mixing171 avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk) 172 avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk) 168 DO jk = 2, jpkm1 !* update momentum & tracer diffusivity with tidal mixing 169 DO jj = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 170 DO ji = 1, jpi 171 avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk) 172 avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk) 173 173 END DO 174 174 END DO 175 175 END DO 176 176 177 DO j j = 2, jpjm1178 DO j i = fs_2, fs_jpim1 ! vector opt.179 DO j k = mikt(ji,jj)+1, jpkm1 !* update momentum & tracer diffusivity with tidal mixing180 avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj ,jk) ) * umask(ji,jj,jk)181 avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji ,jj+1,jk) ) * vmask(ji,jj,jk)177 DO jk = 2, jpkm1 !* update momentum & tracer diffusivity with tidal mixing 178 DO jj = 2, jpjm1 179 DO ji = fs_2, fs_jpim1 ! vector opt. 180 avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj ,jk) ) * wumask(ji,jj,jk) 181 avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji ,jj+1,jk) ) * wvmask(ji,jj,jk) 182 182 END DO 183 183 END DO … … 457 457 ztpc = 0.e0 458 458 zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:) 459 DO j j = 1, jpj460 DO j i = 1, jpi461 DO j k= mikt(ji,jj)+1, jpkm1462 ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)459 DO jk= 2, jpkm1 460 DO jj = 1, jpj 461 DO ji = 1, jpi 462 ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 463 463 END DO 464 464 END DO … … 473 473 zav_tide(:,:,:) = MIN( zav_tide(:,:,:), 60.e-4 ) 474 474 zkz(:,:) = 0.e0 475 DO j j = 1, jpj476 DO j i = 1, jpi477 DO j k = mikt(ji,jj)+1, jpkm1478 zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zav_tide(ji,jj,jk)* tmask(ji,jj,jk)475 DO jk = 2, jpkm1 476 DO jj = 1, jpj 477 DO ji = 1, jpi 478 zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zav_tide(ji,jj,jk)* wmask(ji,jj,jk) 479 479 END DO 480 480 END DO … … 498 498 WRITE(numout,*) ' Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) ) 499 499 500 DO j j = 1, jpj501 DO j i = 1, jpi502 DO j k = mikt(ji,jj)+1, jpkm1503 zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) !kz max = 300 cm2/s500 DO jk = 2, jpkm1 501 DO jj = 1, jpj 502 DO ji = 1, jpi 503 zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk) !kz max = 300 cm2/s 504 504 END DO 505 505 END DO … … 510 510 DO jj = 1, jpj 511 511 DO ji = 1, jpi 512 ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)512 ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 513 513 END DO 514 514 END DO … … 519 519 DO jk = 1, jpk 520 520 ze_z = SUM( e1t(:,:) * e2t(:,:) * zav_tide(:,:,jk) * tmask_i(:,:) ) & 521 & / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )521 & / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) ) 522 522 ztpc = 1.E50 523 523 DO jj = 1, jpj … … 540 540 END DO 541 541 ze_z = SUM( e1t(:,:) * e2t(:,:) * zkz(:,:) * tmask_i(:,:) ) & 542 & / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )542 & / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) ) 543 543 WRITE(numout,*) ' jk= ', jk,' ', ze_z * 1.e4,' cm2/s' 544 544 END DO … … 546 546 zkz(:,:) = az_tmx(:,:,jk) /rn_n2min 547 547 ze_z = SUM( e1t(:,:) * e2t(:,:) * zkz(:,:) * tmask_i(:,:) ) & 548 & / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )548 & / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) ) 549 549 WRITE(numout,*) 550 550 WRITE(numout,*) ' N2 min - jk= ', jk,' ', ze_z * 1.e4,' cm2/s min= ',MINVAL(zkz)*1.e4, & -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/step.F90
r5012 r5098 122 122 IF( lk_zdfkpp ) CALL zdf_kpp( kstp ) ! KPP closure scheme for Kz 123 123 IF( lk_zdfcst ) THEN ! Constant Kz (reset avt, avm[uv] to the background value) 124 avt (:,:,:) = rn_avt0 * tmask(:,:,:)125 avmu(:,:,:) = rn_avm0 * umask(:,:,:)126 avmv(:,:,:) = rn_avm0 * vmask(:,:,:)124 avt (:,:,:) = rn_avt0 * wmask (:,:,:) 125 avmu(:,:,:) = rn_avm0 * wumask(:,:,:) 126 avmv(:,:,:) = rn_avm0 * wvmask(:,:,:) 127 127 ENDIF 128 128 IF( ln_rnf_mouth ) THEN ! increase diffusivity at rivers mouths … … 145 145 ! 146 146 IF( lk_ldfslp ) THEN ! slope of lateral mixing 147 CALL eos( tsb, rhd, gdept_0(:,:,:) ) ! before in situ density 148 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 149 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & ! 150 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 147 CALL eos( tsb, rhd, gdept_0(:,:,:) ) ! before in situ density 148 IF( ln_zps ) CALL zps_hde ( kstp, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 149 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 150 IF( ln_zps .AND. ln_isfcav) & 151 & CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv, & ! Partial steps for top cell (ISF) 152 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 153 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the first ocean level 151 154 IF( ln_traldf_grif ) THEN ! before slope for Griffies operator 152 155 CALL ldf_slp_grif( kstp ) … … 177 180 ! is necessary to compute momentum advection for the rhs of barotropic loop: 178 181 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 179 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient 180 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & ! 181 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 182 IF( ln_zps ) CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient 183 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 184 IF( ln_zps .AND. ln_isfcav) & 185 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps for top cell (ISF) 186 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 187 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 182 188 183 189 ua(:,:,:) = 0.e0 ! set dynamics trends to zero … … 253 259 CALL tra_nxt( kstp ) ! tracer fields at next time step 254 260 CALL eos ( tsa, rhd, rhop, fsdept_n(:,:,:) ) ! Time-filtered in situ density for hpg computation 255 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv, & ! Partial steps: before horizontal gradient 256 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & ! 257 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 261 IF( ln_zps ) CALL zps_hde ( kstp, jpts, tsa, gtsu, gtsv, & ! Partial steps: before horizontal gradient 262 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 263 IF( ln_zps .AND. ln_isfcav) & 264 & CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv, & ! Partial steps for top cell (ISF) 265 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 266 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 258 267 ELSE ! centered hpg (eos then time stepping) 259 268 IF ( .NOT. lk_dynspg_ts ) THEN ! eos already called in time-split case 260 269 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 261 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient 262 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & ! 263 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 270 IF( ln_zps ) CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient 271 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 272 IF( ln_zps .AND. ln_isfcav) & 273 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps for top cell (ISF) 274 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 275 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 264 276 ENDIF 265 277 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90
r4990 r5098 82 82 IF( .NOT. Agrif_Root()) CALL Agrif_Update_Trc( kstp ) ! Update tracer at AGRIF zoom boundaries : children only 83 83 #endif 84 IF( ln_zps ) CALL zps_hde( kstp, jptra, trn, pgtu=gtru, pgtv=gtrv, sgtu=gtrui, sgtv=gtrvi ) ! Partial steps: now horizontal gradient of passive 84 85 IF( ln_zps ) CALL zps_hde ( kstp, jptra, trn, gtru, gtrv ) ! Partial steps: now horizontal gradient of passive 86 IF( ln_zps .AND. ln_isfcav) & 87 & CALL zps_hde_isf( kstp, jptra, trn, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi ) ! Partial steps: now horizontal gradient of passive 85 88 ! tracers at the bottom ocean level 86 89 ! -
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/TOP_SRC/trcini.F90
r4990 r5098 144 144 tra(:,:,:,:) = 0._wp 145 145 IF( ln_zps .AND. .NOT. lk_c1d ) & ! Partial steps: before horizontal gradient of passive 146 & CALL zps_hde( nit000, jptra, trn, pgtu=gtru, pgtv=gtrv, sgtu=gtrui, sgtv=gtrvi ) ! tracers at the bottom ocean level 146 & CALL zps_hde ( nit000, jptra, trn, gtru, gtrv ) ! Partial steps: before horizontal gradient 147 IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) & 148 & CALL zps_hde_isf( nit000, jptra, trn, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi ) ! tracers at the bottom ocean level 149 147 150 148 151 !
Note: See TracChangeset
for help on using the changeset viewer.