Changeset 6006 for branches/2015/dev_MetOffice_merge_2015/NEMOGCM
- Timestamp:
- 2015-12-04T17:56:07+01:00 (9 years ago)
- Location:
- branches/2015/dev_MetOffice_merge_2015/NEMOGCM
- Files:
-
- 23 edited
- 3 copied
-
CONFIG/SHARED/field_def.xml (modified) (1 diff)
-
CONFIG/SHARED/namelist_ref (modified) (4 diffs)
-
NEMO/OFF_SRC/domrea.F90 (modified) (2 diffs)
-
NEMO/OPA_SRC/DIA/diafwb.F90 (modified) (1 diff)
-
NEMO/OPA_SRC/DIA/diaharm.F90 (modified) (4 diffs)
-
NEMO/OPA_SRC/DIA/diahsb.F90 (modified) (10 diffs)
-
NEMO/OPA_SRC/DOM/dom_oce.F90 (modified) (4 diffs)
-
NEMO/OPA_SRC/DOM/domain.F90 (modified) (4 diffs)
-
NEMO/OPA_SRC/DOM/dommsk.F90 (modified) (4 diffs)
-
NEMO/OPA_SRC/DOM/domngb.F90 (modified) (3 diffs)
-
NEMO/OPA_SRC/DOM/domvvl.F90 (modified) (5 diffs)
-
NEMO/OPA_SRC/DOM/domzgr.F90 (modified) (12 diffs)
-
NEMO/OPA_SRC/DOM/iscplhsb.F90 (copied) (copied from branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90)
-
NEMO/OPA_SRC/DOM/iscplini.F90 (copied) (copied from branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90)
-
NEMO/OPA_SRC/DOM/iscplrst.F90 (copied) (copied from branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90)
-
NEMO/OPA_SRC/DOM/istate.F90 (modified) (2 diffs)
-
NEMO/OPA_SRC/DYN/divhor.F90 (modified) (2 diffs)
-
NEMO/OPA_SRC/DYN/dynnxt.F90 (modified) (2 diffs)
-
NEMO/OPA_SRC/IOM/restart.F90 (modified) (2 diffs)
-
NEMO/OPA_SRC/LBC/lbclnk.F90 (modified) (5 diffs)
-
NEMO/OPA_SRC/LBC/lib_mpp.F90 (modified) (2 diffs)
-
NEMO/OPA_SRC/LBC/mppini_2.h90 (modified) (1 diff)
-
NEMO/OPA_SRC/SBC/sbcisf.F90 (modified) (7 diffs)
-
NEMO/OPA_SRC/TRA/tradmp.F90 (modified) (1 diff)
-
NEMO/OPA_SRC/TRA/trasbc.F90 (modified) (6 diffs)
-
NEMO/OPA_SRC/lib_fortran.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/SHARED/field_def.xml
r6005 r6006 200 200 201 201 <!-- * variable related to ice shelf forcing * --> 202 <field id="fwfisf" long_name="Ice shelf melting" unit="kg/m2/s" /> 203 <field id="qisf" long_name="Ice Shelf Heat Flux" unit="W/m2" /> 204 <field id="isfgammat" long_name="transfert coefficient for isf (temperature)" unit="m/s" /> 205 <field id="isfgammas" long_name="transfert coefficient for isf (salinity)" unit="m/s" /> 206 <field id="stbl" long_name="salinity in the Losh tbl" unit="1e-3" /> 207 <field id="ttbl" long_name="temperature in the Losh tbl" unit="degC" /> 202 <field id="fwfisf" long_name="Ice shelf melting" unit="Kg/m2/s" /> 203 <field id="qisf" long_name="Ice Shelf Heat Flux" unit="W/m2" /> 204 <field id="isfgammat" long_name="transfert coefficient for isf (temperature) " unit="m/s" /> 205 <field id="isfgammas" long_name="transfert coefficient for isf (salinity) " unit="m/s" /> 206 <field id="stbl" long_name="salinity in the Losh tbl " unit="PSU" /> 207 <field id="ttbl" long_name="temperature in the Losh tbl " unit="C" /> 208 <field id="utbl" long_name="zonal current in the Losh tbl at T point " unit="m/s" /> 209 <field id="vtbl" long_name="merid current in the Losh tbl at T point " unit="m/s" /> 210 <field id="thermald" long_name="thermal driving of ice shelf melting " unit="C" /> 211 <field id="tfrz" long_name="top freezing point (used to compute melt) " unit="C" /> 212 <field id="tinsitu" long_name="top insitu temperature (used to cmpt melt) " unit="C" /> 213 <field id="ustar" long_name="ustar at T point used in ice shelf melting " unit="m/s" /> 208 214 209 215 <!-- *_oce variables available with ln_blk_clio or ln_blk_core --> -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/SHARED/namelist_ref
r6005 r6006 43 43 cn_ocerst_out = "restart" ! suffix of ocean restart name (output) 44 44 cn_ocerst_outdir = "." ! directory in which to write output ocean restarts 45 ln_iscpl = .false. ! cavity evolution forcing or coupling to ice sheet model 45 46 nn_istate = 0 ! output the initial state (1) or not (0) 46 47 ln_rst_list = .false. ! output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F) … … 128 129 nn_msh = 1 ! create (=1) a mesh file or not (=0) 129 130 rn_hmin = -3. ! min depth of the ocean (>0) or min number of ocean level (<0) 131 rn_isfhmin = 1.00 ! treshold (m) to discriminate grounding ice to floating ice 130 132 rn_e3zps_min= 20. ! partial step thickness is set larger than the minimum of 131 133 rn_e3zps_rat= 0.1 ! rn_e3zps_min and rn_e3zps_rat*e3t, with 0<rn_e3zps_rat<1 … … 212 214 !! namsbc_rnf river runoffs 213 215 !! namsbc_isf ice shelf melting/freezing 216 !! namsbc_iscpl coupling option between land ice model and ocean 214 217 !! namsbc_apr Atmospheric Pressure 215 218 !! namsbc_ssr sea surface restoring term (for T and/or S) … … 456 459 / 457 460 !----------------------------------------------------------------------- 461 &namsbc_iscpl ! land ice / ocean coupling option 462 !----------------------------------------------------------------------- 463 nn_drown = 10 ! number of iteration of the extrapolation loop (fill the new wet cells) 464 ln_hsb = .false. ! activate conservation module (conservation exact after a time of rn_fiscpl) 465 nn_fiscpl = 43800 ! (number of time step) conservation period (maybe should be fix to the coupling frequencey of restart frequency) 466 / 467 !----------------------------------------------------------------------- 458 468 &namsbc_apr ! Atmospheric pressure used as ocean forcing or in bulk 459 469 !----------------------------------------------------------------------- -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OFF_SRC/domrea.F90
r5836 r6006 116 116 & nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl, & 117 117 & nn_it000, nn_itend , nn_date0 , nn_leapy , nn_istate , nn_stock , & 118 & nn_write, ln_ dimgnnn, ln_mskland , ln_cfmeta , ln_clobber, nn_chunksz, nn_euler119 NAMELIST/namdom/ nn_bathy , rn_bathy, rn_e3zps_min, rn_e3zps_rat, nn_msh , rn_hmin, &118 & nn_write, ln_iscpl, ln_dimgnnn, ln_mskland , ln_cfmeta , ln_clobber, nn_chunksz, nn_euler 119 NAMELIST/namdom/ nn_bathy , rn_bathy, rn_e3zps_min, rn_e3zps_rat, nn_msh , rn_hmin, rn_isfhmin, & 120 120 & nn_acc , rn_atfp , rn_rdt , rn_rdtmin , & 121 121 & rn_rdtmax, rn_rdth , nn_baro , nn_closea , ln_crs, & … … 813 813 DO jj = 1, jpjm1 814 814 DO ji = 1, fs_jpim1 ! vector loop 815 umask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:)))816 vmask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:)))815 ssumask(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:))) 816 ssvmask(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:))) 817 817 END DO 818 818 DO ji = 1, jpim1 ! NO vector opt. 819 fmask_i(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) &819 ssfmask(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) & 820 820 & * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 821 821 END DO 822 822 END DO 823 CALL lbc_lnk( umask_i, 'U', 1._wp ) ! Lateral boundary conditions824 CALL lbc_lnk( vmask_i, 'V', 1._wp )825 CALL lbc_lnk( fmask_i, 'F', 1._wp )823 CALL lbc_lnk( ssumask, 'U', 1._wp ) ! Lateral boundary conditions 824 CALL lbc_lnk( ssvmask, 'V', 1._wp ) 825 CALL lbc_lnk( ssfmask, 'F', 1._wp ) 826 826 827 827 ! 3. Ocean/land mask at wu-, wv- and w points -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90
r5836 r6006 460 460 ENDIF 461 461 462 IF( nn_timing == 1 ) CALL timing_st art('dia_fwb')462 IF( nn_timing == 1 ) CALL timing_stop('dia_fwb') 463 463 464 464 9005 FORMAT(1X,A,ES24.16) -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
r5930 r6006 135 135 DO jk=1,nb_ana 136 136 DO ji=1,jpmax_harmo 137 IF (TRIM(tname(jk)) .eq.Wave(ji)%cname_tide) THEN137 IF (TRIM(tname(jk)) == Wave(ji)%cname_tide) THEN 138 138 name(jk) = ji 139 139 EXIT … … 195 195 ! Elevation 196 196 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*tmask_i(ji,jj) 197 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)* umask_i(ji,jj)198 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)* vmask_i(ji,jj)197 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj) 198 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj) 199 199 END DO 200 200 END DO … … 324 324 X1= ana_amp(ji,jj,jh,1) 325 325 X2=-ana_amp(ji,jj,jh,2) 326 out_u(ji,jj, jh) = X1 * umask_i(ji,jj)327 out_u(ji,jj,nb_ana+jh) = X2 * umask_i(ji,jj)326 out_u(ji,jj, jh) = X1 * ssumask(ji,jj) 327 out_u(ji,jj,nb_ana+jh) = X2 * ssumask(ji,jj) 328 328 ENDDO 329 329 ENDDO … … 488 488 DO jj_sd = ji_sd, ninco 489 489 zval2 = ABS(ztmp3(ji_sd,jj_sd)) 490 IF( zval2 .GE.zval1 )THEN490 IF( zval2 >= zval1 )THEN 491 491 ipivot(ji_sd) = jj_sd 492 492 zval1 = zval2 -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
r5643 r6006 46 46 REAL(wp) :: frc_wn_t, frc_wn_s ! global forcing trends 47 47 ! 48 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: surf , ssh_ini ! 48 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: surf 49 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: surf_ini , ssh_ini ! 49 50 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: ssh_hc_loc_ini, ssh_sc_loc_ini ! 50 51 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: hc_loc_ini, sc_loc_ini, e3t_ini ! … … 138 139 ! 2 - Content variations ! 139 140 ! ------------------------ ! 141 ! glob_sum_full is needed because you keep the full interior domain to compute the sum (iscpl) 140 142 zdiff_v2 = 0._wp 141 143 zdiff_hc = 0._wp … … 143 145 144 146 ! volume variation (calculated with ssh) 145 zdiff_v1 = glob_sum ( surf(:,:) * ( sshn(:,:) - ssh_ini(:,:)) )147 zdiff_v1 = glob_sum_full( surf(:,:) * sshn(:,:) - surf_ini(:,:) * ssh_ini(:,:) ) 146 148 147 149 ! heat & salt content variation (associated with ssh) … … 158 160 z2d1(:,:) = surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) 159 161 END IF 160 z_ssh_hc = glob_sum ( z2d0 )161 z_ssh_sc = glob_sum ( z2d1 )162 z_ssh_hc = glob_sum_full( z2d0 ) 163 z_ssh_sc = glob_sum_full( z2d1 ) 162 164 ENDIF 163 165 164 166 DO jk = 1, jpkm1 165 167 ! volume variation (calculated with scale factors) 166 zdiff_v2 = zdiff_v2 + glob_sum ( surf(:,:) * tmask(:,:,jk) &167 & * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk)) )168 ! heat content variation 169 zdiff_hc = zdiff_hc + glob_sum ( surf(:,:) * tmask(:,:,jk) &170 & * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - hc_loc_ini(:,:,jk)) )168 zdiff_v2 = zdiff_v2 + glob_sum_full( surf (:,:) * fse3t_n (:,:,jk) * tmask(:,:,jk) & 169 & - surf_ini(:,:) * e3t_ini(:,:,jk) ) 170 ! heat content variation 171 zdiff_hc = zdiff_hc + glob_sum_full( surf (:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) & 172 & - surf_ini(:,:) * hc_loc_ini(:,:,jk) ) 171 173 ! salt content variation 172 zdiff_sc = zdiff_sc + glob_sum ( surf(:,:) * tmask(:,:,jk) &173 & * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - sc_loc_ini(:,:,jk)) )174 zdiff_sc = zdiff_sc + glob_sum_full( surf (:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) & 175 & - surf_ini(:,:) * sc_loc_ini(:,:,jk) ) 174 176 ENDDO 175 177 … … 191 193 zvol_tot = 0._wp ! total ocean volume (calculated with scale factors) 192 194 DO jk = 1, jpkm1 193 zvol_tot = zvol_tot + glob_sum ( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) )195 zvol_tot = zvol_tot + glob_sum_full( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) ) 194 196 END DO 195 197 … … 214 216 CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot) ! Salt content variation (psu) 215 217 CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp ) ! Heat content variation (1.e20 J) 216 CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9 )! Salt content variation (psu*km3)217 CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9) ! volume ssh variation (km3)218 CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9 ) ! Salt content variation (psu*km3) 219 CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9 ) ! volume ssh variation (km3) 218 220 CALL iom_put( 'bgfrcvol' , frc_v * 1.e-9 ) ! vol - surface forcing (km3) 219 221 CALL iom_put( 'bgfrctem' , frc_t / zvol_tot ) ! hc - surface forcing (C) … … 261 263 CALL iom_get( numror, 'frc_wn_s', frc_wn_s ) 262 264 ENDIF 265 CALL iom_get( numror, jpdom_autoglo, 'surf_ini', surf_ini ) ! ice sheet coupling 263 266 CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini ) 264 267 CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini ) … … 273 276 IF(lwp) WRITE(numout,*) ' dia_hsb at initial state ' 274 277 IF(lwp) WRITE(numout,*) '~~~~~~~' 275 ssh_ini(:,:) = sshn(:,:) ! initial ssh 278 surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:) ! initial ocean surface 279 ssh_ini(:,:) = sshn(:,:) ! initial ssh 276 280 DO jk = 1, jpk 277 e3t_ini (:,:,jk) = fse3t_n(:,:,jk) ! initial vertical scale factors 278 hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk) ! initial heat content 279 sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk) ! initial salt content 281 ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance). 282 e3t_ini (:,:,jk) = fse3t_n(:,:,jk) * tmask(:,:,jk) ! initial vertical scale factors 283 hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk) * tmask(:,:,jk) ! initial heat content 284 sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk) * tmask(:,:,jk) ! initial salt content 280 285 END DO 281 286 frc_v = 0._wp ! volume trend due to forcing … … 312 317 CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s ) 313 318 ENDIF 319 CALL iom_rstput( kt, nitrst, numrow, 'surf_ini', surf_ini ) ! ice sheet coupling 314 320 CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini', ssh_ini ) 315 321 CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini', e3t_ini ) … … 379 385 ! 1 - Allocate memory ! 380 386 ! ------------------- ! 381 ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), &382 & e3t_ini(jpi,jpj,jpk), surf(jpi,jpj), ssh_ini(jpi,jpj), STAT=ierror )387 ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), surf_ini(jpi,jpj), & 388 & e3t_ini(jpi,jpj,jpk), surf(jpi,jpj), ssh_ini(jpi,jpj), STAT=ierror ) 383 389 IF( ierror > 0 ) THEN 384 390 CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' ) ; RETURN -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r5997 r6006 32 32 REAL(wp), PUBLIC :: rn_bathy !: depth of flat bottom (active if nn_bathy=0; if =0 depth=jpkm1) 33 33 REAL(wp), PUBLIC :: rn_hmin !: minimum ocean depth (>0) or minimum number of ocean levels (<0) 34 REAL(wp), PUBLIC :: rn_isfhmin !: threshold to discriminate grounded ice to floating ice 34 35 REAL(wp), PUBLIC :: rn_e3zps_min !: miminum thickness for partial steps (meters) 35 36 REAL(wp), PUBLIC :: rn_e3zps_rat !: minimum thickness ration for partial steps … … 43 44 INTEGER , PUBLIC :: nn_closea !: =0 suppress closed sea/lake from the ORCA domain or not (=1) 44 45 INTEGER , PUBLIC :: nn_euler !: =0 start with forward time step or not (=1) 46 LOGICAL , PUBLIC :: ln_iscpl !: coupling with ice sheet 45 47 LOGICAL , PUBLIC :: ln_crs !: Apply grid coarsening to dynamical model output or online passive tracers 46 48 … … 251 253 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt !: vertical index of the bottom last T- ocean level 252 254 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbku, mbkv !: vertical index of the bottom last U- and W- ocean level 253 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters) 254 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i, umask_i, vmask_i, fmask_i !: interior domain T-point mask 255 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bmask !: land/ocean mask of barotropic stream function 255 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters) 256 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i !: interior domain T-point mask 257 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_h !: internal domain T-point mask (Figure 8.5 NEMO book) 258 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bmask !: land/ocean mask of barotropic stream function 256 259 257 260 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfdep !: top first ocean level (ISF) 258 261 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mikt, miku, mikv, mikf !: first wet T-, U-, V-, F- ocean level (ISF) 259 262 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfdep !: Iceshelf draft (ISF) 260 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask !: surface domain T-point mask 261 263 264 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask, ssumask, ssvmask, ssfmask !: surface mask at T-,U-, V- and F-pts 262 265 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask !: land/ocean mask at T-, U-, V- and F-pts 263 266 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask !: land/ocean mask at WT-, WU- and WV-pts … … 386 389 & hift (jpi,jpj) , hifu (jpi,jpj) , rx1(jpi,jpj) , STAT=ierr(8) ) 387 390 388 ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , & 389 & tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), & 390 & bmask (jpi,jpj) , & 391 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) ) 391 ALLOCATE( mbathy(jpi,jpj) , bathy (jpi,jpj) , & 392 & tmask_i(jpi,jpj) , tmask_h(jpi, jpj), & 393 & ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , & 394 & bmask(jpi,jpj) , & 395 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj) , STAT=ierr(9) ) 392 396 393 397 ! (ISF) Allocation of basic array 394 ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj), 395 & 396 & mikf(jpi,jpj), ssmask(jpi,jpj), STAT=ierr(10) )398 ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj), & 399 & mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) , & 400 & mikf(jpi,jpj), STAT=ierr(10) ) 397 401 398 402 ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk), & -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r5997 r6006 112 112 END DO 113 113 ! ! Inverse of the local depth 114 hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)115 hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)114 hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:) 115 hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:) 116 116 117 117 CALL dom_stp ! time step 118 IF( nmsh /= 0 ) CALL dom_wri ! Create a domain file 118 ! 119 IF( nmsh /= 0 .AND. .NOT. ln_iscpl ) CALL dom_wri ! Create a domain file 120 IF( nmsh /= 0 .AND. ln_iscpl .AND. .NOT. ln_rstart ) CALL dom_wri ! Create a domain file 121 ! 119 122 IF( .NOT.ln_rstart ) CALL dom_ctl ! Domain control 120 123 ! … … 135 138 !!---------------------------------------------------------------------- 136 139 USE ioipsl 137 NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list, &138 nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl ,&139 & nn_it000, nn_itend , nn_date0 , nn_time0 , nn_leapy , nn_istate ,&140 & nn_stock, nn_write , ln_dimgnnn , ln_mskland , ln_clobber, nn_chunksz,&141 & nn_euler, ln_cfmeta 142 NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, 143 & nn_acc , rn_atfp , rn_rdt , rn_rdtmin ,&144 & rn_rdtmax, rn_rdth , nn_closea , ln_crs,&145 & jphgr_msh, &146 & ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, &147 & ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, &140 NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list, & 141 nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl , & 142 & nn_it000, nn_itend , nn_date0 , nn_time0 , nn_leapy , nn_istate , & 143 & nn_stock, nn_write , ln_dimgnnn , ln_mskland , ln_clobber, nn_chunksz, & 144 & nn_euler, ln_cfmeta, ln_iscpl 145 NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, rn_isfhmin, & 146 & nn_acc , rn_atfp , rn_rdt , rn_rdtmin , & 147 & rn_rdtmax, rn_rdth , nn_closea , ln_crs , & 148 & jphgr_msh, & 149 & ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, & 150 & ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, & 148 151 & ppa2, ppkth2, ppacr2 149 152 #if defined key_netcdf4 … … 193 196 WRITE(numout,*) ' overwrite an existing file ln_clobber = ', ln_clobber 194 197 WRITE(numout,*) ' NetCDF chunksize (bytes) nn_chunksz = ', nn_chunksz 198 WRITE(numout,*) ' IS coupling at the restart step ln_iscpl = ', ln_iscpl 195 199 ENDIF 196 200 … … 260 264 WRITE(numout,*) ' min depth of the ocean (>0) or rn_hmin = ', rn_hmin 261 265 WRITE(numout,*) ' min number of ocean level (<0) ' 266 WRITE(numout,*) ' treshold to open the isf cavity rn_isfhmin = ', rn_isfhmin, ' (m)' 262 267 WRITE(numout,*) ' minimum thickness of partial rn_e3zps_min = ', rn_e3zps_min, ' (m)' 263 268 WRITE(numout,*) ' step level rn_e3zps_rat = ', rn_e3zps_rat -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r5930 r6006 183 183 ! -------------------- 184 184 tmask_i(:,:) = ssmask(:,:) ! (ISH) tmask_i = 1 even on the ice shelf 185 186 tmask_h(:,:) = 1._wp ! 0 on the halo and 1 elsewhere 185 187 iif = jpreci ! ??? 186 188 iil = nlci - jpreci + 1 … … 188 190 ijl = nlcj - jprecj + 1 189 191 190 tmask_ i( 1 :iif, : ) = 0._wp ! first columns191 tmask_ i(iil:jpi, : ) = 0._wp ! last columns (including mpp extra columns)192 tmask_ i( : , 1 :ijf) = 0._wp ! first rows193 tmask_ i( : ,ijl:jpj) = 0._wp ! last rows (including mpp extra rows)192 tmask_h( 1 :iif, : ) = 0._wp ! first columns 193 tmask_h(iil:jpi, : ) = 0._wp ! last columns (including mpp extra columns) 194 tmask_h( : , 1 :ijf) = 0._wp ! first rows 195 tmask_h( : ,ijl:jpj) = 0._wp ! last rows (including mpp extra rows) 194 196 195 197 ! north fold mask … … 202 204 IF( mjg(nlej) == jpjglo ) THEN ! only half of the nlcj-1 row 203 205 DO ji = iif+1, iil-1 204 tmask_ i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))206 tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji)) 205 207 END DO 206 208 ENDIF 207 209 ENDIF 210 211 tmask_i(:,:) = tmask_i(:,:) * tmask_h(:,:) 212 208 213 IF( jperio == 5 .OR. jperio == 6 ) THEN ! F-point pivot 209 214 tpol( 1 :jpiglo) = 0._wp … … 225 230 END DO 226 231 END DO 227 ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point232 ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet cell at u point 228 233 DO jj = 1, jpjm1 229 234 DO ji = 1, fs_jpim1 ! vector loop 230 umask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:)))231 vmask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:)))235 ssumask(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:))) 236 ssvmask(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:))) 232 237 END DO 233 238 DO ji = 1, jpim1 ! NO vector opt. 234 fmask_i(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) &239 ssfmask(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) & 235 240 & * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 236 241 END DO 237 242 END DO 238 CALL lbc_lnk( umask , 'U', 1._wp ) ! Lateral boundary conditions239 CALL lbc_lnk( vmask , 'V', 1._wp )240 CALL lbc_lnk( fmask , 'F', 1._wp )241 CALL lbc_lnk( umask_i, 'U', 1._wp ) ! Lateral boundary conditions242 CALL lbc_lnk( vmask_i, 'V', 1._wp )243 CALL lbc_lnk( fmask_i, 'F', 1._wp )243 CALL lbc_lnk( umask , 'U', 1._wp ) ! Lateral boundary conditions 244 CALL lbc_lnk( vmask , 'V', 1._wp ) 245 CALL lbc_lnk( fmask , 'F', 1._wp ) 246 CALL lbc_lnk( ssumask, 'U', 1._wp ) ! Lateral boundary conditions 247 CALL lbc_lnk( ssvmask, 'V', 1._wp ) 248 CALL lbc_lnk( ssfmask, 'F', 1._wp ) 244 249 245 250 ! 3. Ocean/land mask at wu-, wv- and w points -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90
r5836 r6006 28 28 CONTAINS 29 29 30 SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid )30 SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid, kkk ) 31 31 !!---------------------------------------------------------------------- 32 32 !! *** ROUTINE dom_ngb *** … … 39 39 REAL(wp) , INTENT(in ) :: plon, plat ! longitude,latitude of the point 40 40 INTEGER , INTENT( out) :: kii, kjj ! i-,j-index of the closes grid point 41 INTEGER , INTENT(in ), OPTIONAL :: kkk ! k-index of the mask level used 41 42 CHARACTER(len=1), INTENT(in ) :: cdgrid ! grid name 'T', 'U', 'V', 'W' 42 43 ! 44 INTEGER :: ik ! working level 43 45 INTEGER , DIMENSION(2) :: iloc 44 46 REAL(wp) :: zlon, zmini … … 51 53 ! 52 54 zmask(:,:) = 0._wp 55 ik = 1 56 IF ( PRESENT(kkk) ) ik=kkk 53 57 SELECT CASE( cdgrid ) 54 CASE( 'U' ) ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej, 1)55 CASE( 'V' ) ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej, 1)56 CASE( 'F' ) ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej, 1)57 CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej, 1)58 CASE( 'U' ) ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,ik) 59 CASE( 'V' ) ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,ik) 60 CASE( 'F' ) ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,ik) 61 CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,ik) 58 62 END SELECT 59 63 60 zlon = MOD( plon + 720., 360. ) ! plon between 0 and 360 61 zglam(:,:) = MOD( zglam(:,:) + 720., 360. ) ! glam between 0 and 360 62 IF( zlon > 270. ) zlon = zlon - 360. ! zlon between -90 and 270 63 IF( zlon < 90. ) WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360. ! glam between -180 and 180 64 IF (jphgr_msh /= 2 .AND. jphgr_msh /= 3) THEN 65 zlon = MOD( plon + 720., 360. ) ! plon between 0 and 360 66 zglam(:,:) = MOD( zglam(:,:) + 720., 360. ) ! glam between 0 and 360 67 IF( zlon > 270. ) zlon = zlon - 360. ! zlon between -90 and 270 68 IF( zlon < 90. ) WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360. ! glam between -180 and 180 69 zglam(:,:) = zglam(:,:) - zlon 70 ELSE 71 zglam(:,:) = zglam(:,:) - plon 72 END IF 64 73 65 zglam(:,:) = zglam(:,:) - zlon66 74 zgphi(:,:) = zgphi(:,:) - plat 67 75 zdist(:,:) = zglam(:,:) * zglam(:,:) + zgphi(:,:) * zgphi(:,:) -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r5836 r6006 24 24 USE in_out_manager ! I/O manager 25 25 USE iom ! I/O manager library 26 USE restart 26 USE restart, ONLY : rst_read_open ! ocean restart 27 27 USE lib_mpp ! distributed memory computing library 28 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 192 192 hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 193 193 END DO 194 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) )195 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) )194 hur_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1. - ssumask(:,:) ) 195 hvr_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1. - ssvmask(:,:) ) 196 196 197 197 ! Restoring frequencies for z_tilde coordinate … … 325 325 ! -------------------------------------------------- 326 326 IF( ln_vvl_ztilde ) THEN 327 IF( kt .GT.nit000 ) THEN327 IF( kt > nit000 ) THEN 328 328 DO jk = 1, jpkm1 329 329 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & … … 419 419 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 420 420 ! - ML - test: for the moment, stop simulation for too large e3_t variations 421 IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT.- rn_zdef_max ) ) THEN421 IF( ( z_tmax > rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN 422 422 IF( lk_mpp ) THEN 423 423 CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) ) … … 538 538 END DO 539 539 ! ! Inverse of the local depth 540 hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)541 hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)540 hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:) 541 hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:) 542 542 543 543 CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5836 r6006 379 379 !! - bathy : meter bathymetry (in meters) 380 380 !!---------------------------------------------------------------------- 381 INTEGER :: ji, jj, j l, jk ! dummy loop indices381 INTEGER :: ji, jj, jk ! dummy loop indices 382 382 INTEGER :: inum ! temporary logical unit 383 383 INTEGER :: ierror ! error flag … … 541 541 CALL iom_close( inum ) 542 542 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 543 544 ! set grounded point to 0 545 ! (a treshold could be set here if needed, or set it offline based on the grounded fraction) 546 WHERE ( bathy(:,:) <= risfdep(:,:) + rn_isfhmin ) 547 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 548 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 549 END WHERE 543 550 END IF 544 551 ! … … 578 585 ! 579 586 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 580 ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin581 IF ( ln_isfcav ) THEN582 WHERE (bathy == risfdep)583 bathy = 0.0_wp ; risfdep = 0.0_wp584 END WHERE585 END IF586 ! end patch587 587 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level 588 588 ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth … … 835 835 END SUBROUTINE zgr_bot_level 836 836 837 838 !!---------------------------------------------------------------------- 839 !! *** ROUTINE zgr_ bot_level ***837 SUBROUTINE zgr_top_level 838 !!---------------------------------------------------------------------- 839 !! *** ROUTINE zgr_top_level *** 840 840 !! 841 841 !! ** Purpose : defines the vertical index of ocean top (mik. arrays) … … 963 963 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 964 964 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 965 REAL(wp) :: zmax ! Maximum depth966 965 REAL(wp) :: zdiff ! temporary scalar 967 REAL(wp) :: z refdep! temporary scalar966 REAL(wp) :: zmax ! temporary scalar 968 967 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 969 968 !!--------------------------------------------------------------------- … … 1004 1003 END DO 1005 1004 1006 IF ( ln_isfcav ) CALL zgr_isf1007 1008 1005 ! Scale factors and depth at T- and W-points 1009 1006 DO jk = 1, jpk ! intitialization to the reference z-coordinate … … 1013 1010 e3w_0 (:,:,jk) = e3w_1d (jk) 1014 1011 END DO 1012 1013 ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf 1014 IF ( ln_isfcav ) CALL zgr_isf 1015 1016 ! Scale factors and depth at T- and W-points 1017 IF ( .NOT. ln_isfcav ) THEN 1018 DO jj = 1, jpj 1019 DO ji = 1, jpi 1020 ik = mbathy(ji,jj) 1021 IF( ik > 0 ) THEN ! ocean point only 1022 ! max ocean level case 1023 IF( ik == jpkm1 ) THEN 1024 zdepwp = bathy(ji,jj) 1025 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1026 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1027 e3t_0(ji,jj,ik ) = ze3tp 1028 e3t_0(ji,jj,ik+1) = ze3tp 1029 e3w_0(ji,jj,ik ) = ze3wp 1030 e3w_0(ji,jj,ik+1) = ze3tp 1031 gdepw_0(ji,jj,ik+1) = zdepwp 1032 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1033 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1034 ! 1035 ELSE ! standard case 1036 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1037 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1038 ENDIF 1039 !gm Bug? check the gdepw_1d 1040 ! ... on ik 1041 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1042 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1043 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1044 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1045 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1046 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1047 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1048 ! ... on ik+1 1049 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1050 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1051 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 1052 ENDIF 1053 ENDIF 1054 END DO 1055 END DO 1056 ! 1057 it = 0 1058 DO jj = 1, jpj 1059 DO ji = 1, jpi 1060 ik = mbathy(ji,jj) 1061 IF( ik > 0 ) THEN ! ocean point only 1062 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1063 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1064 ! test 1065 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1066 IF( zdiff <= 0._wp .AND. lwp ) THEN 1067 it = it + 1 1068 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1069 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1070 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1071 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1072 ENDIF 1073 ENDIF 1074 END DO 1075 END DO 1076 END IF 1077 ! 1078 ! Scale factors and depth at U-, V-, UW and VW-points 1079 DO jk = 1, jpk ! initialisation to z-scale factors 1080 e3u_0 (:,:,jk) = e3t_1d(jk) 1081 e3v_0 (:,:,jk) = e3t_1d(jk) 1082 e3uw_0(:,:,jk) = e3w_1d(jk) 1083 e3vw_0(:,:,jk) = e3w_1d(jk) 1084 END DO 1085 1086 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 1087 DO jj = 1, jpjm1 1088 DO ji = 1, fs_jpim1 ! vector opt. 1089 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 1090 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 1091 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 1092 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 1093 END DO 1094 END DO 1095 END DO 1096 IF ( ln_isfcav ) THEN 1097 ! (ISF) define e3uw (adapted for 2 cells in the water column) 1098 DO jj = 2, jpjm1 1099 DO ji = 2, fs_jpim1 ! vector opt. 1100 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 1101 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 1102 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) & 1103 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) ) 1104 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 1105 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 1106 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) & 1107 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) ) 1108 END DO 1109 END DO 1110 END IF 1111 1112 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 1113 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 1114 ! 1115 1116 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1117 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk) 1118 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk) 1119 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk) 1120 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk) 1121 END DO 1122 1123 ! Scale factor at F-point 1124 DO jk = 1, jpk ! initialisation to z-scale factors 1125 e3f_0(:,:,jk) = e3t_1d(jk) 1126 END DO 1127 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors 1128 DO jj = 1, jpjm1 1129 DO ji = 1, fs_jpim1 ! vector opt. 1130 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 1131 END DO 1132 END DO 1133 END DO 1134 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions 1135 ! 1136 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1137 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk) 1138 END DO 1139 !!gm bug ? : must be a do loop with mj0,mj1 1015 1140 ! 1141 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 1142 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 1143 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 1144 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 1145 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 1146 1147 ! Control of the sign 1148 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' ) 1149 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' ) 1150 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' ) 1151 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' ) 1152 1153 ! Compute gdep3w_0 (vertical sum of e3w) 1154 IF ( ln_isfcav ) THEN ! if cavity 1155 WHERE (misfdep == 0) misfdep = 1 1156 DO jj = 1,jpj 1157 DO ji = 1,jpi 1158 gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1159 DO jk = 2, misfdep(ji,jj) 1160 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1161 END DO 1162 IF (misfdep(ji,jj) >= 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 1163 DO jk = misfdep(ji,jj) + 1, jpk 1164 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1165 END DO 1166 END DO 1167 END DO 1168 ELSE ! no cavity 1169 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1170 DO jk = 2, jpk 1171 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1172 END DO 1173 END IF 1174 ! ! ================= ! 1175 IF(lwp .AND. ll_print) THEN ! Control print ! 1176 ! ! ================= ! 1177 DO jj = 1,jpj 1178 DO ji = 1, jpi 1179 ik = MAX( mbathy(ji,jj), 1 ) 1180 zprt(ji,jj,1) = e3t_0 (ji,jj,ik) 1181 zprt(ji,jj,2) = e3w_0 (ji,jj,ik) 1182 zprt(ji,jj,3) = e3u_0 (ji,jj,ik) 1183 zprt(ji,jj,4) = e3v_0 (ji,jj,ik) 1184 zprt(ji,jj,5) = e3f_0 (ji,jj,ik) 1185 zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 1186 END DO 1187 END DO 1188 WRITE(numout,*) 1189 WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1190 WRITE(numout,*) 1191 WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1192 WRITE(numout,*) 1193 WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1194 WRITE(numout,*) 1195 WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1196 WRITE(numout,*) 1197 WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1198 WRITE(numout,*) 1199 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1200 ENDIF 1201 ! 1202 CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 1203 ! 1204 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps') 1205 ! 1206 END SUBROUTINE zgr_zps 1207 1208 SUBROUTINE zgr_isf 1209 !!---------------------------------------------------------------------- 1210 !! *** ROUTINE zgr_isf *** 1211 !! 1212 !! ** Purpose : check the bathymetry in levels 1213 !! 1214 !! ** Method : THe water column have to contained at least 2 cells 1215 !! Bathymetry and isfdraft are modified (dig/close) to respect 1216 !! this criterion. 1217 !! 1218 !! ** Action : - test compatibility between isfdraft and bathy 1219 !! - bathy and isfdraft are modified 1220 !!---------------------------------------------------------------------- 1221 !! 1222 INTEGER :: ji, jj, jl, jk ! dummy loop indices 1223 INTEGER :: ik, it ! temporary integers 1224 INTEGER :: icompt, ibtest ! (ISF) 1225 INTEGER :: ibtestim1, ibtestip1 ! (ISF) 1226 INTEGER :: ibtestjm1, ibtestjp1 ! (ISF) 1227 REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t 1228 REAL(wp) :: zmax ! Maximum and minimum depth 1229 REAL(wp) :: zbathydiff ! isf temporary scalar 1230 REAL(wp) :: zrisfdepdiff ! isf temporary scalar 1231 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 1232 REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t 1233 REAL(wp) :: zdiff ! temporary scalar 1234 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH) 1235 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH) 1236 !!--------------------------------------------------------------------- 1237 ! 1238 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1239 ! 1240 CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 1241 CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy ) 1242 1243 ! (ISF) compute misfdep 1244 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ; misfdep(:,:) = 1 ! open water : set misfdep to 1 1245 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level 1246 END WHERE 1247 1248 ! Compute misfdep for ocean points (i.e. first wet level) 1249 ! find the first ocean level such that the first level thickness 1250 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 1251 ! e3t_0 is the reference level thickness 1252 DO jk = 2, jpkm1 1253 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1254 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1255 END DO 1256 WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) ) 1257 risfdep(:,:) = 0. ; misfdep(:,:) = 1 1258 END WHERE 1259 1260 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation 1261 icompt = 0 1262 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1263 DO jl = 1, 10 1264 ! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments) 1265 WHERE (bathy(:,:) <= risfdep(:,:) + rn_isfhmin) 1266 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1267 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 1268 END WHERE 1269 WHERE (mbathy(:,:) <= 0) 1270 misfdep(:,:) = 0; risfdep(:,:) = 0._wp 1271 mbathy (:,:) = 0; bathy (:,:) = 0._wp 1272 ENDWHERE 1273 IF( lk_mpp ) THEN 1274 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1275 CALL lbc_lnk( zbathy, 'T', 1. ) 1276 misfdep(:,:) = INT( zbathy(:,:) ) 1277 1278 CALL lbc_lnk( risfdep,'T', 1. ) 1279 CALL lbc_lnk( bathy, 'T', 1. ) 1280 1281 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1282 CALL lbc_lnk( zbathy, 'T', 1. ) 1283 mbathy(:,:) = INT( zbathy(:,:) ) 1284 ENDIF 1285 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1286 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west 1287 misfdep(jpi,:) = misfdep( 2 ,:) 1288 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 1289 mbathy(jpi,:) = mbathy( 2 ,:) 1290 ENDIF 1291 1292 ! split last cell if possible (only where water column is 2 cell or less) 1293 ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss). 1294 IF ( .NOT. ln_iscpl) THEN 1295 DO jk = jpkm1, 1, -1 1296 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1297 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 1298 mbathy(:,:) = jk 1299 bathy(:,:) = zmax 1300 END WHERE 1301 END DO 1302 END IF 1303 1304 ! split top cell if possible (only where water column is 2 cell or less) 1305 DO jk = 2, jpkm1 1306 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1307 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 1308 misfdep(:,:) = jk 1309 risfdep(:,:) = zmax 1310 END WHERE 1311 END DO 1312 1313 1314 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 1315 DO jj = 1, jpj 1316 DO ji = 1, jpi 1317 ! find the minimum change option: 1318 ! test bathy 1319 IF (risfdep(ji,jj) > 1) THEN 1320 IF ( .NOT. ln_iscpl ) THEN 1321 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1322 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1323 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1324 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1325 1326 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1327 IF (zbathydiff <= zrisfdepdiff) THEN 1328 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 1329 mbathy(ji,jj)= mbathy(ji,jj) + 1 1330 ELSE 1331 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1332 misfdep(ji,jj) = misfdep(ji,jj) - 1 1333 END IF 1334 ENDIF 1335 ELSE 1336 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1337 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1338 misfdep(ji,jj) = misfdep(ji,jj) - 1 1339 END IF 1340 END IF 1341 END IF 1342 END DO 1343 END DO 1344 1345 ! At least 2 levels for water thickness at T, U, and V point. 1346 DO jj = 1, jpj 1347 DO ji = 1, jpi 1348 ! find the minimum change option: 1349 ! test bathy 1350 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1351 IF ( .NOT. ln_iscpl ) THEN 1352 zbathydiff =ABS(bathy(ji,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 1353 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1354 zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj) ) & 1355 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1356 IF (zbathydiff <= zrisfdepdiff) THEN 1357 mbathy(ji,jj) = mbathy(ji,jj) + 1 1358 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1359 ELSE 1360 misfdep(ji,jj)= misfdep(ji,jj) - 1 1361 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1362 END IF 1363 ELSE 1364 misfdep(ji,jj)= misfdep(ji,jj) - 1 1365 risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1366 END IF 1367 ENDIF 1368 END DO 1369 END DO 1370 1371 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1) 1372 DO jj = 1, jpjm1 1373 DO ji = 1, jpim1 1374 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1375 IF ( .NOT. ln_iscpl ) THEN 1376 zbathydiff =ABS(bathy(ji,jj ) - ( gdepw_1d(mbathy (ji,jj)+1) & 1377 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat ))) 1378 zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) & 1379 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 1380 IF (zbathydiff <= zrisfdepdiff) THEN 1381 mbathy(ji,jj) = mbathy(ji,jj) + 1 1382 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat ) 1383 ELSE 1384 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1385 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1386 END IF 1387 ELSE 1388 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1389 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1390 END IF 1391 ENDIF 1392 END DO 1393 END DO 1394 1395 IF( lk_mpp ) THEN 1396 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1397 CALL lbc_lnk( zbathy, 'T', 1. ) 1398 misfdep(:,:) = INT( zbathy(:,:) ) 1399 1400 CALL lbc_lnk( risfdep,'T', 1. ) 1401 CALL lbc_lnk( bathy, 'T', 1. ) 1402 1403 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1404 CALL lbc_lnk( zbathy, 'T', 1. ) 1405 mbathy(:,:) = INT( zbathy(:,:) ) 1406 ENDIF 1407 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1) 1408 DO jj = 1, jpjm1 1409 DO ji = 1, jpim1 1410 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN 1411 IF ( .NOT. ln_iscpl ) THEN 1412 zbathydiff =ABS( bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) & 1413 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 1414 zrisfdepdiff=ABS(risfdep(ji,jj ) - ( gdepw_1d(misfdep(ji,jj ) ) & 1415 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat ))) 1416 IF (zbathydiff <= zrisfdepdiff) THEN 1417 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 1418 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 1419 ELSE 1420 misfdep(ji,jj) = misfdep(ji,jj) - 1 1421 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1422 END IF 1423 ELSE 1424 misfdep(ji,jj) = misfdep(ji,jj) - 1 1425 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1426 END IF 1427 ENDIF 1428 END DO 1429 END DO 1430 1431 1432 IF( lk_mpp ) THEN 1433 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1434 CALL lbc_lnk( zbathy, 'T', 1. ) 1435 misfdep(:,:) = INT( zbathy(:,:) ) 1436 1437 CALL lbc_lnk( risfdep,'T', 1. ) 1438 CALL lbc_lnk( bathy, 'T', 1. ) 1439 1440 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1441 CALL lbc_lnk( zbathy, 'T', 1. ) 1442 mbathy(:,:) = INT( zbathy(:,:) ) 1443 ENDIF 1444 1445 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1) 1446 DO jj = 1, jpjm1 1447 DO ji = 1, jpim1 1448 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1449 IF ( .NOT. ln_iscpl ) THEN 1450 zbathydiff =ABS( bathy(ji ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 1451 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat ))) 1452 zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) & 1453 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 1454 IF (zbathydiff <= zrisfdepdiff) THEN 1455 mbathy(ji,jj) = mbathy(ji,jj) + 1 1456 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1457 ELSE 1458 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1459 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1460 END IF 1461 ELSE 1462 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1463 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1464 ENDIF 1465 ENDIF 1466 ENDDO 1467 ENDDO 1468 1469 IF( lk_mpp ) THEN 1470 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1471 CALL lbc_lnk( zbathy, 'T', 1. ) 1472 misfdep(:,:) = INT( zbathy(:,:) ) 1473 1474 CALL lbc_lnk( risfdep,'T', 1. ) 1475 CALL lbc_lnk( bathy, 'T', 1. ) 1476 1477 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1478 CALL lbc_lnk( zbathy, 'T', 1. ) 1479 mbathy(:,:) = INT( zbathy(:,:) ) 1480 ENDIF 1481 1482 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1) 1483 DO jj = 1, jpjm1 1484 DO ji = 1, jpim1 1485 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN 1486 IF ( .NOT. ln_iscpl ) THEN 1487 zbathydiff =ABS( bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) & 1488 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 1489 zrisfdepdiff=ABS(risfdep(ji ,jj) - ( gdepw_1d(misfdep(ji ,jj) ) & 1490 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat ))) 1491 IF (zbathydiff <= zrisfdepdiff) THEN 1492 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1 1493 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 1494 ELSE 1495 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1496 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1497 END IF 1498 ELSE 1499 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1500 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1501 ENDIF 1502 ENDIF 1503 ENDDO 1504 ENDDO 1505 1506 IF( lk_mpp ) THEN 1507 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1508 CALL lbc_lnk( zbathy, 'T', 1. ) 1509 misfdep(:,:) = INT( zbathy(:,:) ) 1510 1511 CALL lbc_lnk( risfdep,'T', 1. ) 1512 CALL lbc_lnk( bathy, 'T', 1. ) 1513 1514 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1515 CALL lbc_lnk( zbathy, 'T', 1. ) 1516 mbathy(:,:) = INT( zbathy(:,:) ) 1517 ENDIF 1518 END DO 1519 ! end dig bathy/ice shelf to be compatible 1520 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 1521 DO jl = 1,20 1522 1523 ! remove single point "bay" on isf coast line in the ice shelf draft' 1524 DO jk = 2, jpk 1525 WHERE (misfdep==0) misfdep=jpk 1526 zmask=0._wp 1527 WHERE (misfdep <= jk) zmask=1 1528 DO jj = 2, jpjm1 1529 DO ji = 2, jpim1 1530 IF (misfdep(ji,jj) == jk) THEN 1531 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1532 IF (ibtest <= 1) THEN 1533 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 1534 IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk 1535 END IF 1536 END IF 1537 END DO 1538 END DO 1539 END DO 1540 WHERE (misfdep==jpk) 1541 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1542 END WHERE 1543 IF( lk_mpp ) THEN 1544 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1545 CALL lbc_lnk( zbathy, 'T', 1. ) 1546 misfdep(:,:) = INT( zbathy(:,:) ) 1547 1548 CALL lbc_lnk( risfdep,'T', 1. ) 1549 CALL lbc_lnk( bathy, 'T', 1. ) 1550 1551 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1552 CALL lbc_lnk( zbathy, 'T', 1. ) 1553 mbathy(:,:) = INT( zbathy(:,:) ) 1554 ENDIF 1555 1556 ! remove single point "bay" on bathy coast line beneath an ice shelf' 1557 DO jk = jpk,1,-1 1558 zmask=0._wp 1559 WHERE (mbathy >= jk ) zmask=1 1560 DO jj = 2, jpjm1 1561 DO ji = 2, jpim1 1562 IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN 1563 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1564 IF (ibtest <= 1) THEN 1565 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 1566 IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0 1567 END IF 1568 END IF 1569 END DO 1570 END DO 1571 END DO 1572 WHERE (mbathy==0) 1573 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1574 END WHERE 1575 IF( lk_mpp ) THEN 1576 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1577 CALL lbc_lnk( zbathy, 'T', 1. ) 1578 misfdep(:,:) = INT( zbathy(:,:) ) 1579 1580 CALL lbc_lnk( risfdep,'T', 1. ) 1581 CALL lbc_lnk( bathy, 'T', 1. ) 1582 1583 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1584 CALL lbc_lnk( zbathy, 'T', 1. ) 1585 mbathy(:,:) = INT( zbathy(:,:) ) 1586 ENDIF 1587 1588 ! fill hole in ice shelf 1589 zmisfdep = misfdep 1590 zrisfdep = risfdep 1591 WHERE (zmisfdep <= 1._wp) zmisfdep=jpk 1592 DO jj = 2, jpjm1 1593 DO ji = 2, jpim1 1594 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj ) 1595 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1) 1596 IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj ) ) ibtestim1 = jpk 1597 IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj ) ) ibtestip1 = jpk 1598 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj-1) ) ibtestjm1 = jpk 1599 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj+1) ) ibtestjp1 = jpk 1600 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1601 IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN 1602 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 1603 END IF 1604 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN 1605 misfdep(ji,jj) = ibtest 1606 risfdep(ji,jj) = gdepw_1d(ibtest) 1607 ENDIF 1608 ENDDO 1609 ENDDO 1610 1611 IF( lk_mpp ) THEN 1612 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1613 CALL lbc_lnk( zbathy, 'T', 1. ) 1614 misfdep(:,:) = INT( zbathy(:,:) ) 1615 1616 CALL lbc_lnk( risfdep, 'T', 1. ) 1617 CALL lbc_lnk( bathy, 'T', 1. ) 1618 1619 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1620 CALL lbc_lnk( zbathy, 'T', 1. ) 1621 mbathy(:,:) = INT( zbathy(:,:) ) 1622 ENDIF 1623 ! 1624 !! fill hole in bathymetry 1625 zmbathy (:,:)=mbathy (:,:) 1626 DO jj = 2, jpjm1 1627 DO ji = 2, jpim1 1628 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj ) 1629 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1) 1630 IF( zmbathy(ji,jj) < misfdep(ji-1,jj ) ) ibtestim1 = 0 1631 IF( zmbathy(ji,jj) < misfdep(ji+1,jj ) ) ibtestip1 = 0 1632 IF( zmbathy(ji,jj) < misfdep(ji ,jj-1) ) ibtestjm1 = 0 1633 IF( zmbathy(ji,jj) < misfdep(ji ,jj+1) ) ibtestjp1 = 0 1634 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1635 IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN 1636 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 1637 END IF 1638 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN 1639 mbathy(ji,jj) = ibtest 1640 bathy(ji,jj) = gdepw_1d(ibtest+1) 1641 ENDIF 1642 END DO 1643 END DO 1644 IF( lk_mpp ) THEN 1645 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1646 CALL lbc_lnk( zbathy, 'T', 1. ) 1647 misfdep(:,:) = INT( zbathy(:,:) ) 1648 1649 CALL lbc_lnk( risfdep, 'T', 1. ) 1650 CALL lbc_lnk( bathy, 'T', 1. ) 1651 1652 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1653 CALL lbc_lnk( zbathy, 'T', 1. ) 1654 mbathy(:,:) = INT( zbathy(:,:) ) 1655 ENDIF 1656 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1657 DO jj = 1, jpjm1 1658 DO ji = 1, jpim1 1659 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1660 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1661 END IF 1662 END DO 1663 END DO 1664 IF( lk_mpp ) THEN 1665 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1666 CALL lbc_lnk( zbathy, 'T', 1. ) 1667 misfdep(:,:) = INT( zbathy(:,:) ) 1668 1669 CALL lbc_lnk( risfdep, 'T', 1. ) 1670 CALL lbc_lnk( bathy, 'T', 1. ) 1671 1672 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1673 CALL lbc_lnk( zbathy, 'T', 1. ) 1674 mbathy(:,:) = INT( zbathy(:,:) ) 1675 ENDIF 1676 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1677 DO jj = 1, jpjm1 1678 DO ji = 1, jpim1 1679 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1680 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ; 1681 END IF 1682 END DO 1683 END DO 1684 IF( lk_mpp ) THEN 1685 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1686 CALL lbc_lnk( zbathy, 'T', 1. ) 1687 misfdep(:,:) = INT( zbathy(:,:) ) 1688 1689 CALL lbc_lnk( risfdep,'T', 1. ) 1690 CALL lbc_lnk( bathy, 'T', 1. ) 1691 1692 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1693 CALL lbc_lnk( zbathy, 'T', 1. ) 1694 mbathy(:,:) = INT( zbathy(:,:) ) 1695 ENDIF 1696 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1697 DO jj = 1, jpjm1 1698 DO ji = 1, jpi 1699 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1700 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1701 END IF 1702 END DO 1703 END DO 1704 IF( lk_mpp ) THEN 1705 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1706 CALL lbc_lnk( zbathy, 'T', 1. ) 1707 misfdep(:,:) = INT( zbathy(:,:) ) 1708 1709 CALL lbc_lnk( risfdep,'T', 1. ) 1710 CALL lbc_lnk( bathy, 'T', 1. ) 1711 1712 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1713 CALL lbc_lnk( zbathy, 'T', 1. ) 1714 mbathy(:,:) = INT( zbathy(:,:) ) 1715 ENDIF 1716 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1717 DO jj = 1, jpjm1 1718 DO ji = 1, jpi 1719 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1720 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ; 1721 END IF 1722 END DO 1723 END DO 1724 IF( lk_mpp ) THEN 1725 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1726 CALL lbc_lnk( zbathy, 'T', 1. ) 1727 misfdep(:,:) = INT( zbathy(:,:) ) 1728 1729 CALL lbc_lnk( risfdep,'T', 1. ) 1730 CALL lbc_lnk( bathy, 'T', 1. ) 1731 1732 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1733 CALL lbc_lnk( zbathy, 'T', 1. ) 1734 mbathy(:,:) = INT( zbathy(:,:) ) 1735 ENDIF 1736 ! if not compatible after all check, mask T 1737 DO jj = 1, jpj 1738 DO ji = 1, jpi 1739 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 1740 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ; 1741 END IF 1742 END DO 1743 END DO 1744 1745 WHERE (mbathy(:,:) == 1) 1746 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 1747 END WHERE 1748 END DO 1749 ! end check compatibility ice shelf/bathy 1750 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1751 IF( icompt == 0 ) THEN 1752 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry' 1753 ELSE 1754 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 1755 ENDIF 1756 1757 ! compute scale factor and depth at T- and W- points 1016 1758 DO jj = 1, jpj 1017 1759 DO ji = 1, jpi … … 1035 1777 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1036 1778 ENDIF 1779 ! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1037 1780 !gm Bug? check the gdepw_1d 1038 1781 ! ... on ik … … 1040 1783 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1041 1784 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1042 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1043 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1044 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1045 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1785 e3t_0 (ji,jj,ik ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik ) 1786 e3w_0 (ji,jj,ik ) = gdept_0(ji,jj,ik ) - gdept_1d(ik-1) 1046 1787 ! ... on ik+1 1047 1788 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1048 1789 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1049 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)1050 1790 ENDIF 1051 1791 ENDIF … … 1073 1813 END DO 1074 1814 ! 1075 IF ( ln_isfcav ) THEN1076 1815 ! (ISF) Definition of e3t, u, v, w for ISF case 1077 1078 1079 1080 1081 1082 1816 DO jj = 1, jpj 1817 DO ji = 1, jpi 1818 ik = misfdep(ji,jj) 1819 IF( ik > 1 ) THEN ! ice shelf point only 1820 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1821 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1083 1822 !gm Bug? check the gdepw_0 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 e3w_0 (ji,jj,ik ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))1096 1823 ! ... on ik 1824 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1825 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1826 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1827 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1828 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1829 1830 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1831 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1832 ENDIF 1833 ! ... on ik / ik-1 1834 e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik) 1835 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1097 1836 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1098 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1837 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1838 ENDIF 1839 END DO 1840 END DO 1841 1842 it = 0 1843 DO jj = 1, jpj 1844 DO ji = 1, jpi 1845 ik = misfdep(ji,jj) 1846 IF( ik > 1 ) THEN ! ice shelf point only 1847 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1848 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1849 ! test 1850 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1851 IF( zdiff <= 0. .AND. lwp ) THEN 1852 it = it + 1 1853 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1854 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1855 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1856 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1099 1857 ENDIF 1100 END DO1858 ENDIF 1101 1859 END DO 1102 !1103 it = 01104 DO jj = 1, jpj1105 DO ji = 1, jpi1106 ik = misfdep(ji,jj)1107 IF( ik > 1 ) THEN ! ice shelf point only1108 e3tp (ji,jj) = e3t_0(ji,jj,ik )1109 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )1110 ! test1111 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik )1112 IF( zdiff <= 0. .AND. lwp ) THEN1113 it = it + 11114 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj1115 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)1116 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff1117 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj)1118 ENDIF1119 ENDIF1120 END DO1121 END DO1122 END IF1123 ! END (ISF)1124 1125 ! Scale factors and depth at U-, V-, UW and VW-points1126 DO jk = 1, jpk ! initialisation to z-scale factors1127 e3u_0 (:,:,jk) = e3t_1d(jk)1128 e3v_0 (:,:,jk) = e3t_1d(jk)1129 e3uw_0(:,:,jk) = e3w_1d(jk)1130 e3vw_0(:,:,jk) = e3w_1d(jk)1131 END DO1132 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors1133 DO jj = 1, jpjm11134 DO ji = 1, fs_jpim1 ! vector opt.1135 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )1136 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )1137 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )1138 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )1139 END DO1140 END DO1141 END DO1142 IF ( ln_isfcav ) THEN1143 ! (ISF) define e3uw (adapted for 2 cells in the water column)1144 DO jj = 2, jpjm11145 DO ji = 2, fs_jpim1 ! vector opt.1146 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))1147 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))1148 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) &1149 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) )1150 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))1151 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))1152 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) &1153 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) )1154 END DO1155 END DO1156 END IF1157 1158 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions1159 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp )1160 !1161 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1162 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)1163 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)1164 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)1165 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)1166 END DO1167 1168 ! Scale factor at F-point1169 DO jk = 1, jpk ! initialisation to z-scale factors1170 e3f_0(:,:,jk) = e3t_1d(jk)1171 END DO1172 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors1173 DO jj = 1, jpjm11174 DO ji = 1, fs_jpim1 ! vector opt.1175 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )1176 END DO1177 END DO1178 END DO1179 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions1180 !1181 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1182 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)1183 END DO1184 !!gm bug ? : must be a do loop with mj0,mj11185 !1186 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 21187 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)1188 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)1189 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)1190 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)1191 1192 ! Control of the sign1193 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )1194 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )1195 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )1196 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )1197 1198 ! Compute gdep3w_0 (vertical sum of e3w)1199 IF ( ln_isfcav ) THEN ! if cavity1200 WHERE (misfdep == 0) misfdep = 11201 DO jj = 1,jpj1202 DO ji = 1,jpi1203 gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)1204 DO jk = 2, misfdep(ji,jj)1205 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1206 END DO1207 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))1208 DO jk = misfdep(ji,jj) + 1, jpk1209 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1210 END DO1211 END DO1212 END DO1213 ELSE ! no cavity1214 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)1215 DO jk = 2, jpk1216 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)1217 END DO1218 END IF1219 ! ! ================= !1220 IF(lwp .AND. ll_print) THEN ! Control print !1221 ! ! ================= !1222 DO jj = 1,jpj1223 DO ji = 1, jpi1224 ik = MAX( mbathy(ji,jj), 1 )1225 zprt(ji,jj,1) = e3t_0 (ji,jj,ik)1226 zprt(ji,jj,2) = e3w_0 (ji,jj,ik)1227 zprt(ji,jj,3) = e3u_0 (ji,jj,ik)1228 zprt(ji,jj,4) = e3v_0 (ji,jj,ik)1229 zprt(ji,jj,5) = e3f_0 (ji,jj,ik)1230 zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)1231 END DO1232 END DO1233 WRITE(numout,*)1234 WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1235 WRITE(numout,*)1236 WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1237 WRITE(numout,*)1238 WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1239 WRITE(numout,*)1240 WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1241 WRITE(numout,*)1242 WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1243 WRITE(numout,*)1244 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1245 ENDIF1246 !1247 CALL wrk_dealloc( jpi, jpj, jpk, zprt )1248 !1249 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps')1250 !1251 END SUBROUTINE zgr_zps1252 1253 SUBROUTINE zgr_isf1254 !!----------------------------------------------------------------------1255 !! *** ROUTINE zgr_isf ***1256 !!1257 !! ** Purpose : check the bathymetry in levels1258 !!1259 !! ** Method : THe water column have to contained at least 2 cells1260 !! Bathymetry and isfdraft are modified (dig/close) to respect1261 !! this criterion.1262 !!1263 !!1264 !! ** Action : - test compatibility between isfdraft and bathy1265 !! - bathy and isfdraft are modified1266 !!----------------------------------------------------------------------1267 !!1268 INTEGER :: ji, jj, jk, jl ! dummy loop indices1269 INTEGER :: ik, it ! temporary integers1270 INTEGER :: id, jd, nprocd1271 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF)1272 LOGICAL :: ll_print ! Allow control print for debugging1273 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points1274 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t1275 REAL(wp) :: zmax, zmin ! Maximum and minimum depth1276 REAL(wp) :: zdiff ! temporary scalar1277 REAL(wp) :: zrefdep ! temporary scalar1278 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar1279 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH)1280 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH)1281 !!---------------------------------------------------------------------1282 !1283 IF( nn_timing == 1 ) CALL timing_start('zgr_isf')1284 !1285 CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)1286 CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )1287 1288 1289 ! (ISF) compute misfdep1290 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ; misfdep(:,:) = 1 ! open water : set misfdep to 11291 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level1292 END WHERE1293 1294 ! Compute misfdep for ocean points (i.e. first wet level)1295 ! find the first ocean level such that the first level thickness1296 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where1297 ! e3t_0 is the reference level thickness1298 DO jk = 2, jpkm11299 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1300 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+11301 1860 END DO 1302 WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp)1303 risfdep(:,:) = 0. ; misfdep(:,:) = 11304 END WHERE1305 1306 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation1307 icompt = 01308 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together1309 DO jl = 1, 101310 WHERE (bathy(:,:) .EQ. risfdep(:,:) )1311 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp1312 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp1313 END WHERE1314 WHERE (mbathy(:,:) <= 0)1315 misfdep(:,:) = 0; risfdep(:,:) = 0._wp1316 mbathy (:,:) = 0; bathy (:,:) = 0._wp1317 ENDWHERE1318 IF( lk_mpp ) THEN1319 zbathy(:,:) = FLOAT( misfdep(:,:) )1320 CALL lbc_lnk( zbathy, 'T', 1. )1321 misfdep(:,:) = INT( zbathy(:,:) )1322 CALL lbc_lnk( risfdep, 'T', 1. )1323 CALL lbc_lnk( bathy, 'T', 1. )1324 zbathy(:,:) = FLOAT( mbathy(:,:) )1325 CALL lbc_lnk( zbathy, 'T', 1. )1326 mbathy(:,:) = INT( zbathy(:,:) )1327 ENDIF1328 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN1329 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west1330 misfdep(jpi,:) = misfdep( 2 ,:)1331 ENDIF1332 1333 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN1334 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west1335 mbathy(jpi,:) = mbathy( 2 ,:)1336 ENDIF1337 1338 ! split last cell if possible (only where water column is 2 cell or less)1339 DO jk = jpkm1, 1, -11340 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1341 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)1342 mbathy(:,:) = jk1343 bathy(:,:) = zmax1344 END WHERE1345 END DO1346 1347 ! split top cell if possible (only where water column is 2 cell or less)1348 DO jk = 2, jpkm11349 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1350 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)1351 misfdep(:,:) = jk1352 risfdep(:,:) = zmax1353 END WHERE1354 END DO1355 1356 1357 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition1358 DO jj = 1, jpj1359 DO ji = 1, jpi1360 ! find the minimum change option:1361 ! test bathy1362 IF (risfdep(ji,jj) .GT. 1) THEN1363 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) &1364 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))1365 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &1366 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))1367 1368 IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN1369 IF (zbathydiff .LE. zrisfdepdiff) THEN1370 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )1371 mbathy(ji,jj)= mbathy(ji,jj) + 11372 ELSE1373 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )1374 misfdep(ji,jj) = misfdep(ji,jj) - 11375 END IF1376 END IF1377 END IF1378 END DO1379 END DO1380 1381 ! At least 2 levels for water thickness at T, U, and V point.1382 DO jj = 1, jpj1383 DO ji = 1, jpi1384 ! find the minimum change option:1385 ! test bathy1386 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1387 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1)&1388 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))1389 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &1390 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))1391 IF (zbathydiff .LE. zrisfdepdiff) THEN1392 mbathy(ji,jj) = mbathy(ji,jj) + 11393 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )1394 ELSE1395 misfdep(ji,jj)= misfdep(ji,jj) - 11396 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )1397 END IF1398 ENDIF1399 END DO1400 END DO1401 1402 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)1403 DO jj = 1, jpjm11404 DO ji = 1, jpim11405 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1406 zbathydiff =ABS(bathy(ji,jj ) - (gdepw_1d(mbathy (ji,jj)+1) &1407 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat )))1408 zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) &1409 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))1410 IF (zbathydiff .LE. zrisfdepdiff) THEN1411 mbathy(ji,jj) = mbathy(ji,jj) + 11412 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) &1413 & + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat )1414 ELSE1415 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 11416 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) &1417 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )1418 END IF1419 ENDIF1420 END DO1421 END DO1422 1423 IF( lk_mpp ) THEN1424 zbathy(:,:) = FLOAT( misfdep(:,:) )1425 CALL lbc_lnk( zbathy, 'T', 1. )1426 misfdep(:,:) = INT( zbathy(:,:) )1427 CALL lbc_lnk( risfdep, 'T', 1. )1428 CALL lbc_lnk( bathy, 'T', 1. )1429 zbathy(:,:) = FLOAT( mbathy(:,:) )1430 CALL lbc_lnk( zbathy, 'T', 1. )1431 mbathy(:,:) = INT( zbathy(:,:) )1432 ENDIF1433 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)1434 DO jj = 1, jpjm11435 DO ji = 1, jpim11436 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN1437 zbathydiff =ABS( bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) &1438 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))1439 zrisfdepdiff=ABS(risfdep(ji,jj ) - (gdepw_1d(misfdep(ji,jj ) ) &1440 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat )))1441 IF (zbathydiff .LE. zrisfdepdiff) THEN1442 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 11443 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )1444 ELSE1445 misfdep(ji,jj) = misfdep(ji,jj) - 11446 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat )1447 END IF1448 ENDIF1449 END DO1450 END DO1451 1452 1453 IF( lk_mpp ) THEN1454 zbathy(:,:) = FLOAT( misfdep(:,:) )1455 CALL lbc_lnk( zbathy, 'T', 1. )1456 misfdep(:,:) = INT( zbathy(:,:) )1457 CALL lbc_lnk( risfdep, 'T', 1. )1458 CALL lbc_lnk( bathy, 'T', 1. )1459 zbathy(:,:) = FLOAT( mbathy(:,:) )1460 CALL lbc_lnk( zbathy, 'T', 1. )1461 mbathy(:,:) = INT( zbathy(:,:) )1462 ENDIF1463 1464 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)1465 DO jj = 1, jpjm11466 DO ji = 1, jpim11467 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1468 zbathydiff =ABS( bathy(ji ,jj) - (gdepw_1d(mbathy (ji,jj)+1) &1469 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat )))1470 zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) &1471 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))1472 IF (zbathydiff .LE. zrisfdepdiff) THEN1473 mbathy(ji,jj) = mbathy(ji,jj) + 11474 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )1475 ELSE1476 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 11477 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )1478 END IF1479 ENDIF1480 ENDDO1481 ENDDO1482 1483 IF( lk_mpp ) THEN1484 zbathy(:,:) = FLOAT( misfdep(:,:) )1485 CALL lbc_lnk( zbathy, 'T', 1. )1486 misfdep(:,:) = INT( zbathy(:,:) )1487 CALL lbc_lnk( risfdep, 'T', 1. )1488 CALL lbc_lnk( bathy, 'T', 1. )1489 zbathy(:,:) = FLOAT( mbathy(:,:) )1490 CALL lbc_lnk( zbathy, 'T', 1. )1491 mbathy(:,:) = INT( zbathy(:,:) )1492 ENDIF1493 1494 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1)1495 DO jj = 1, jpjm11496 DO ji = 1, jpim11497 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1498 zbathydiff =ABS( bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) &1499 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))1500 zrisfdepdiff=ABS(risfdep(ji ,jj) - (gdepw_1d(misfdep(ji ,jj) ) &1501 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat )))1502 IF (zbathydiff .LE. zrisfdepdiff) THEN1503 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 11504 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) &1505 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )1506 ELSE1507 misfdep(ji,jj) = misfdep(ji ,jj) - 11508 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) &1509 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat )1510 END IF1511 ENDIF1512 ENDDO1513 ENDDO1514 1515 IF( lk_mpp ) THEN1516 zbathy(:,:) = FLOAT( misfdep(:,:) )1517 CALL lbc_lnk( zbathy, 'T', 1. )1518 misfdep(:,:) = INT( zbathy(:,:) )1519 CALL lbc_lnk( risfdep, 'T', 1. )1520 CALL lbc_lnk( bathy, 'T', 1. )1521 zbathy(:,:) = FLOAT( mbathy(:,:) )1522 CALL lbc_lnk( zbathy, 'T', 1. )1523 mbathy(:,:) = INT( zbathy(:,:) )1524 ENDIF1525 END DO1526 ! end dig bathy/ice shelf to be compatible1527 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness1528 DO jl = 1,201529 1530 ! remove single point "bay" on isf coast line in the ice shelf draft'1531 DO jk = 2, jpk1532 WHERE (misfdep==0) misfdep=jpk1533 zmask=01534 WHERE (misfdep .LE. jk) zmask=11535 DO jj = 2, jpjm11536 DO ji = 2, jpim11537 IF (misfdep(ji,jj) .EQ. jk) THEN1538 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)1539 IF (ibtest .LE. 1) THEN1540 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+11541 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk1542 END IF1543 END IF1544 END DO1545 END DO1546 END DO1547 WHERE (misfdep==jpk)1548 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.1549 END WHERE1550 IF( lk_mpp ) THEN1551 zbathy(:,:) = FLOAT( misfdep(:,:) )1552 CALL lbc_lnk( zbathy, 'T', 1. )1553 misfdep(:,:) = INT( zbathy(:,:) )1554 CALL lbc_lnk( risfdep, 'T', 1. )1555 CALL lbc_lnk( bathy, 'T', 1. )1556 zbathy(:,:) = FLOAT( mbathy(:,:) )1557 CALL lbc_lnk( zbathy, 'T', 1. )1558 mbathy(:,:) = INT( zbathy(:,:) )1559 ENDIF1560 1561 ! remove single point "bay" on bathy coast line beneath an ice shelf'1562 DO jk = jpk,1,-11563 zmask=01564 WHERE (mbathy .GE. jk ) zmask=11565 DO jj = 2, jpjm11566 DO ji = 2, jpim11567 IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN1568 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)1569 IF (ibtest .LE. 1) THEN1570 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-11571 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 01572 END IF1573 END IF1574 END DO1575 END DO1576 END DO1577 WHERE (mbathy==0)1578 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.1579 END WHERE1580 IF( lk_mpp ) THEN1581 zbathy(:,:) = FLOAT( misfdep(:,:) )1582 CALL lbc_lnk( zbathy, 'T', 1. )1583 misfdep(:,:) = INT( zbathy(:,:) )1584 CALL lbc_lnk( risfdep, 'T', 1. )1585 CALL lbc_lnk( bathy, 'T', 1. )1586 zbathy(:,:) = FLOAT( mbathy(:,:) )1587 CALL lbc_lnk( zbathy, 'T', 1. )1588 mbathy(:,:) = INT( zbathy(:,:) )1589 ENDIF1590 1591 ! fill hole in ice shelf1592 zmisfdep = misfdep1593 zrisfdep = risfdep1594 WHERE (zmisfdep .LE. 1) zmisfdep=jpk1595 DO jj = 2, jpjm11596 DO ji = 2, jpim11597 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj )1598 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1)1599 IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj ) - 1)1600 IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj ) - 1)1601 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji ,jj-1) - 1)1602 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji ,jj+1) - 1)1603 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)1604 IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN1605 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp1606 END IF1607 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN1608 misfdep(ji,jj) = ibtest1609 risfdep(ji,jj) = gdepw_1d(ibtest)1610 ENDIF1611 ENDDO1612 ENDDO1613 1614 IF( lk_mpp ) THEN1615 zbathy(:,:) = FLOAT( misfdep(:,:) )1616 CALL lbc_lnk( zbathy, 'T', 1. )1617 misfdep(:,:) = INT( zbathy(:,:) )1618 CALL lbc_lnk( risfdep, 'T', 1. )1619 CALL lbc_lnk( bathy, 'T', 1. )1620 zbathy(:,:) = FLOAT( mbathy(:,:) )1621 CALL lbc_lnk( zbathy, 'T', 1. )1622 mbathy(:,:) = INT( zbathy(:,:) )1623 ENDIF1624 !1625 !! fill hole in bathymetry1626 zmbathy (:,:)=mbathy (:,:)1627 DO jj = 2, jpjm11628 DO ji = 2, jpim11629 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj )1630 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1)1631 IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj ) + 1)1632 IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj ) ) ibtestip1 = 01633 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj-1) ) ibtestjm1 = 01634 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj+1) ) ibtestjp1 = 01635 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)1636 IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN1637 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;1638 END IF1639 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN1640 mbathy(ji,jj) = ibtest1641 bathy(ji,jj) = gdepw_1d(ibtest+1)1642 ENDIF1643 END DO1644 END DO1645 IF( lk_mpp ) THEN1646 zbathy(:,:) = FLOAT( misfdep(:,:) )1647 CALL lbc_lnk( zbathy, 'T', 1. )1648 misfdep(:,:) = INT( zbathy(:,:) )1649 CALL lbc_lnk( risfdep, 'T', 1. )1650 CALL lbc_lnk( bathy, 'T', 1. )1651 zbathy(:,:) = FLOAT( mbathy(:,:) )1652 CALL lbc_lnk( zbathy, 'T', 1. )1653 mbathy(:,:) = INT( zbathy(:,:) )1654 ENDIF1655 ! if not compatible after all check (ie U point water column less than 2 cells), mask U1656 DO jj = 1, jpjm11657 DO ji = 1, jpim11658 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN1659 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;1660 END IF1661 END DO1662 END DO1663 IF( lk_mpp ) THEN1664 zbathy(:,:) = FLOAT( misfdep(:,:) )1665 CALL lbc_lnk( zbathy, 'T', 1. )1666 misfdep(:,:) = INT( zbathy(:,:) )1667 CALL lbc_lnk( risfdep, 'T', 1. )1668 CALL lbc_lnk( bathy, 'T', 1. )1669 zbathy(:,:) = FLOAT( mbathy(:,:) )1670 CALL lbc_lnk( zbathy, 'T', 1. )1671 mbathy(:,:) = INT( zbathy(:,:) )1672 ENDIF1673 ! if not compatible after all check (ie U point water column less than 2 cells), mask U1674 DO jj = 1, jpjm11675 DO ji = 1, jpim11676 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN1677 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ;1678 END IF1679 END DO1680 END DO1681 IF( lk_mpp ) THEN1682 zbathy(:,:) = FLOAT( misfdep(:,:) )1683 CALL lbc_lnk( zbathy, 'T', 1. )1684 misfdep(:,:) = INT( zbathy(:,:) )1685 CALL lbc_lnk( risfdep, 'T', 1. )1686 CALL lbc_lnk( bathy, 'T', 1. )1687 zbathy(:,:) = FLOAT( mbathy(:,:) )1688 CALL lbc_lnk( zbathy, 'T', 1. )1689 mbathy(:,:) = INT( zbathy(:,:) )1690 ENDIF1691 ! if not compatible after all check (ie V point water column less than 2 cells), mask V1692 DO jj = 1, jpjm11693 DO ji = 1, jpi1694 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN1695 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;1696 END IF1697 END DO1698 END DO1699 IF( lk_mpp ) THEN1700 zbathy(:,:) = FLOAT( misfdep(:,:) )1701 CALL lbc_lnk( zbathy, 'T', 1. )1702 misfdep(:,:) = INT( zbathy(:,:) )1703 CALL lbc_lnk( risfdep, 'T', 1. )1704 CALL lbc_lnk( bathy, 'T', 1. )1705 zbathy(:,:) = FLOAT( mbathy(:,:) )1706 CALL lbc_lnk( zbathy, 'T', 1. )1707 mbathy(:,:) = INT( zbathy(:,:) )1708 ENDIF1709 ! if not compatible after all check (ie V point water column less than 2 cells), mask V1710 DO jj = 1, jpjm11711 DO ji = 1, jpi1712 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN1713 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ;1714 END IF1715 END DO1716 END DO1717 IF( lk_mpp ) THEN1718 zbathy(:,:) = FLOAT( misfdep(:,:) )1719 CALL lbc_lnk( zbathy, 'T', 1. )1720 misfdep(:,:) = INT( zbathy(:,:) )1721 CALL lbc_lnk( risfdep, 'T', 1. )1722 CALL lbc_lnk( bathy, 'T', 1. )1723 zbathy(:,:) = FLOAT( mbathy(:,:) )1724 CALL lbc_lnk( zbathy, 'T', 1. )1725 mbathy(:,:) = INT( zbathy(:,:) )1726 ENDIF1727 ! if not compatible after all check, mask T1728 DO jj = 1, jpj1729 DO ji = 1, jpi1730 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN1731 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ;1732 END IF1733 END DO1734 END DO1735 1736 WHERE (mbathy(:,:) == 1)1737 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp1738 END WHERE1739 END DO1740 ! end check compatibility ice shelf/bathy1741 ! remove very shallow ice shelf (less than ~ 10m if 75L)1742 WHERE (misfdep(:,:) <= 5)1743 misfdep = 1; risfdep = 0.0_wp;1744 END WHERE1745 1746 IF( icompt == 0 ) THEN1747 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry'1748 ELSE1749 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'1750 ENDIF1751 1861 1752 1862 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) … … 1755 1865 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1756 1866 1757 END SUBROUTINE 1867 END SUBROUTINE zgr_isf 1758 1868 1759 1869 SUBROUTINE zgr_sco … … 2084 2194 fsde3w(:,:,:) = gdep3w_0(:,:,:) 2085 2195 ! 2086 where (e3t_0 (:,:,:) .eq.0.0) e3t_0(:,:,:) = 1.02087 where (e3u_0 (:,:,:) .eq.0.0) e3u_0(:,:,:) = 1.02088 where (e3v_0 (:,:,:) .eq.0.0) e3v_0(:,:,:) = 1.02089 where (e3f_0 (:,:,:) .eq.0.0) e3f_0(:,:,:) = 1.02090 where (e3w_0 (:,:,:) .eq.0.0) e3w_0(:,:,:) = 1.02091 where (e3uw_0 (:,:,:) .eq.0.0) e3uw_0(:,:,:) = 1.02092 where (e3vw_0 (:,:,:) .eq.0.0) e3vw_0(:,:,:) = 1.02196 where (e3t_0 (:,:,:) == 0.0) e3t_0(:,:,:) = 1.0_wp 2197 where (e3u_0 (:,:,:) == 0.0) e3u_0(:,:,:) = 1.0_wp 2198 where (e3v_0 (:,:,:) == 0.0) e3v_0(:,:,:) = 1.0_wp 2199 where (e3f_0 (:,:,:) == 0.0) e3f_0(:,:,:) = 1.0_wp 2200 where (e3w_0 (:,:,:) == 0.0) e3w_0(:,:,:) = 1.0_wp 2201 where (e3uw_0 (:,:,:) == 0.0) e3uw_0(:,:,:) = 1.0_wp 2202 where (e3vw_0 (:,:,:) == 0.0) e3vw_0(:,:,:) = 1.0_wp 2093 2203 2094 2204 #if defined key_agrif -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r5930 r6006 35 35 USE dtauvd ! data: U & V current (dta_uvd routine) 36 36 USE domvvl ! varying vertical mesh 37 USE iscplrst ! ice sheet coupling 37 38 ! 38 39 USE in_out_manager ! I/O manager … … 85 86 IF( ln_rstart ) THEN ! Restart from a file 86 87 ! ! ------------------- 87 CALL rst_read ! Read the restart file 88 CALL day_init ! model calendar (using both namelist and restart infos) 88 CALL rst_read ! Read the restart file 89 IF (ln_iscpl) CALL iscpl_stp ! extraloate restart to wet and dry 90 CALL day_init ! model calendar (using both namelist and restart infos) 89 91 ELSE 90 92 ! ! Start from rest -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90
r5836 r6006 23 23 USE sbcrnf ! river runoff 24 24 USE sbcisf ! ice shelf 25 USE iscplhsb ! ice sheet / ocean coupling 26 USE iscplini ! ice sheet / ocean coupling 25 27 ! 26 28 USE in_out_manager ! I/O manager … … 93 95 IF( ln_divisf .AND. nn_isf > 0 ) CALL sbc_isf_div( hdivn ) !== ice shelf ==! (update hdivn field) 94 96 ! 97 IF( ln_iscpl .AND. ln_hsb ) CALL iscpl_div( hdivn ) !== ice sheet ==! (update hdivn field) 98 ! 95 99 CALL lbc_lnk( hdivn, 'T', 1. ) !== lateral boundary cond. ==! (no sign change) 96 100 ! -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r5930 r6006 95 95 ! 96 96 INTEGER :: ji, jj, jk ! dummy loop indices 97 INTEGER :: iku, ikv ! local integers98 97 REAL(wp) :: zue3a, zue3n, zue3b, zuf ! local scalars 99 98 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - … … 320 319 hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 321 320 END DO 322 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) )323 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) )321 hur_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 322 hvr_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 324 323 ENDIF 325 324 ! -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
r6005 r6006 142 142 CALL iom_rstput( kt, nitrst, numrow, 'sshn' , sshn ) 143 143 CALL iom_rstput( kt, nitrst, numrow, 'rhop' , rhop ) 144 145 ! extra variable needed for the ice sheet coupling 146 IF ( ln_iscpl ) THEN 147 CALL iom_rstput( kt, nitrst, numrow, 'tmask' , tmask ) ! need to extrapolate T/S 148 CALL iom_rstput( kt, nitrst, numrow, 'umask' , umask ) ! need to correct barotropic velocity 149 CALL iom_rstput( kt, nitrst, numrow, 'vmask' , vmask ) ! need to correct barotropic velocity 150 CALL iom_rstput( kt, nitrst, numrow, 'smask' , ssmask ) ! need to correct barotropic velocity 151 CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) ) ! need to compute temperature correction 152 CALL iom_rstput( kt, nitrst, numrow, 'fse3u_n', fse3u_n(:,:,:) ) ! need to compute bt conservation 153 CALL iom_rstput( kt, nitrst, numrow, 'fse3v_n', fse3v_n(:,:,:) ) ! need to compute bt conservation 154 CALL iom_rstput( kt, nitrst, numrow, 'fsdepw_n', fsdepw_n(:,:,:) ) ! need to compute extrapolation if vvl 155 END IF 144 156 ENDIF 145 157 146 158 IF (ln_diurnal) CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst ) 147 159 148 160 IF( kt == nitrst ) THEN 149 161 CALL iom_close( numrow ) ! close the restart file (only at last time step) … … 209 221 REAL(wp) :: zrdt, zrdttra1 210 222 INTEGER :: jk 211 LOGICAL :: llok212 223 !!---------------------------------------------------------------------- 213 224 -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90
r5429 r6006 17 17 !!---------------------------------------------------------------------- 18 18 !! lbc_lnk : generic interface for mpp_lnk_3d and mpp_lnk_2d routines defined in lib_mpp 19 !! lbc_sum : generic interface for mpp_lnk_sum_3d and mpp_lnk_sum_2d routines defined in lib_mpp 19 20 !! lbc_lnk_e : generic interface for mpp_lnk_2d_e routine defined in lib_mpp 20 21 !! lbc_bdy_lnk : generic interface for mpp_lnk_bdy_2d and mpp_lnk_bdy_3d routines defined in lib_mpp … … 31 32 END INTERFACE 32 33 34 !JMM interface not defined if not key_mpp_mpi : likely do not compile without this CPP key !!!! 35 INTERFACE lbc_sum 36 MODULE PROCEDURE mpp_lnk_sum_3d, mpp_lnk_sum_2d 37 END INTERFACE 38 33 39 INTERFACE lbc_bdy_lnk 34 40 MODULE PROCEDURE mpp_lnk_bdy_2d, mpp_lnk_bdy_3d … … 45 51 PUBLIC lbc_lnk ! ocean lateral boundary conditions 46 52 PUBLIC lbc_lnk_multi ! modified ocean lateral boundary conditions 53 PUBLIC lbc_sum 47 54 PUBLIC lbc_lnk_e 48 55 PUBLIC lbc_bdy_lnk ! ocean lateral BDY boundary conditions … … 59 66 !! Default option shared memory computing 60 67 !!---------------------------------------------------------------------- 61 !! lbc_lnk : generic interface for lbc_lnk_3d and lbc_lnk_2d 62 !! lbc_lnk_3d : set the lateral boundary condition on a 3D variable on ocean mesh 63 !! lbc_lnk_2d : set the lateral boundary condition on a 2D variable on ocean mesh 64 !! lbc_bdy_lnk : set the lateral BDY boundary condition 68 !! lbc_sum : generic interface for mpp_lnk_sum_3d and mpp_lnk_sum_2d 69 !! lbc_lnk_sum_3d: compute sum over the halos on a 3D variable on ocean mesh 70 !! lbc_lnk_sum_3d: compute sum over the halos on a 2D variable on ocean mesh 71 !! lbc_lnk : generic interface for lbc_lnk_3d and lbc_lnk_2d 72 !! lbc_lnk_3d : set the lateral boundary condition on a 3D variable on ocean mesh 73 !! lbc_lnk_2d : set the lateral boundary condition on a 2D variable on ocean mesh 74 !! lbc_bdy_lnk : set the lateral BDY boundary condition 65 75 !!---------------------------------------------------------------------- 66 76 USE oce ! ocean dynamics and tracers … … 74 84 INTERFACE lbc_lnk 75 85 MODULE PROCEDURE lbc_lnk_3d_gather, lbc_lnk_3d, lbc_lnk_2d 86 END INTERFACE 87 88 INTERFACE lbc_sum 89 MODULE PROCEDURE mpp_lnk_sum_3d, mpp_lnk_sum_2d 76 90 END INTERFACE 77 91 -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
r5836 r6006 72 72 PUBLIC mpp_lnk_3d, mpp_lnk_3d_gather, mpp_lnk_2d, mpp_lnk_2d_e 73 73 PUBLIC mpp_lnk_2d_9 74 PUBLIC mpp_lnk_sum_3d, mpp_lnk_sum_2d 74 75 PUBLIC mppscatter, mppgather 75 76 PUBLIC mpp_ini_ice, mpp_ini_znl … … 1402 1403 END SUBROUTINE mpp_lnk_2d_e 1403 1404 1405 SUBROUTINE mpp_lnk_sum_3d( ptab, cd_type, psgn, cd_mpp, pval ) 1406 !!---------------------------------------------------------------------- 1407 !! *** routine mpp_lnk_sum_3d *** 1408 !! 1409 !! ** Purpose : Message passing manadgement (sum the overlap region) 1410 !! 1411 !! ** Method : Use mppsend and mpprecv function for passing mask 1412 !! between processors following neighboring subdomains. 1413 !! domain parameters 1414 !! nlci : first dimension of the local subdomain 1415 !! nlcj : second dimension of the local subdomain 1416 !! nbondi : mark for "east-west local boundary" 1417 !! nbondj : mark for "north-south local boundary" 1418 !! noea : number for local neighboring processors 1419 !! nowe : number for local neighboring processors 1420 !! noso : number for local neighboring processors 1421 !! nono : number for local neighboring processors 1422 !! 1423 !! ** Action : ptab with update value at its periphery 1424 !! 1425 !!---------------------------------------------------------------------- 1426 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: ptab ! 3D array on which the boundary condition is applied 1427 CHARACTER(len=1) , INTENT(in ) :: cd_type ! define the nature of ptab array grid-points 1428 ! ! = T , U , V , F , W points 1429 REAL(wp) , INTENT(in ) :: psgn ! =-1 the sign change across the north fold boundary 1430 ! ! = 1. , the sign is kept 1431 CHARACTER(len=3), OPTIONAL , INTENT(in ) :: cd_mpp ! fill the overlap area only 1432 REAL(wp) , OPTIONAL , INTENT(in ) :: pval ! background value (used at closed boundaries) 1433 !! 1434 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1435 INTEGER :: imigr, iihom, ijhom ! temporary integers 1436 INTEGER :: ml_req1, ml_req2, ml_err ! for key_mpi_isend 1437 REAL(wp) :: zland 1438 INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ml_stat ! for key_mpi_isend 1439 ! 1440 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zt3ns, zt3sn ! 3d for north-south & south-north 1441 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zt3ew, zt3we ! 3d for east-west & west-east 1442 1443 !!---------------------------------------------------------------------- 1444 1445 ALLOCATE( zt3ns(jpi,jprecj,jpk,2), zt3sn(jpi,jprecj,jpk,2), & 1446 & zt3ew(jpj,jpreci,jpk,2), zt3we(jpj,jpreci,jpk,2) ) 1447 1448 ! 1449 IF( PRESENT( pval ) ) THEN ; zland = pval ! set land value 1450 ELSE ; zland = 0.e0 ! zero by default 1451 ENDIF 1452 1453 ! 1. standard boundary treatment 1454 ! ------------------------------ 1455 ! 2. East and west directions exchange 1456 ! ------------------------------------ 1457 ! we play with the neigbours AND the row number because of the periodicity 1458 ! 1459 SELECT CASE ( nbondi ) ! Read lateral conditions 1460 CASE ( -1, 0, 1 ) ! all exept 2 (i.e. close case) 1461 iihom = nlci-jpreci 1462 DO jl = 1, jpreci 1463 zt3ew(:,jl,:,1) = ptab(jl ,:,:) ; ptab(jl ,:,:) = 0.0_wp 1464 zt3we(:,jl,:,1) = ptab(iihom+jl,:,:) ; ptab(iihom+jl,:,:) = 0.0_wp 1465 END DO 1466 END SELECT 1467 ! 1468 ! ! Migrations 1469 imigr = jpreci * jpj * jpk 1470 ! 1471 SELECT CASE ( nbondi ) 1472 CASE ( -1 ) 1473 CALL mppsend( 2, zt3we(1,1,1,1), imigr, noea, ml_req1 ) 1474 CALL mpprecv( 1, zt3ew(1,1,1,2), imigr, noea ) 1475 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 1476 CASE ( 0 ) 1477 CALL mppsend( 1, zt3ew(1,1,1,1), imigr, nowe, ml_req1 ) 1478 CALL mppsend( 2, zt3we(1,1,1,1), imigr, noea, ml_req2 ) 1479 CALL mpprecv( 1, zt3ew(1,1,1,2), imigr, noea ) 1480 CALL mpprecv( 2, zt3we(1,1,1,2), imigr, nowe ) 1481 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 1482 IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err) 1483 CASE ( 1 ) 1484 CALL mppsend( 1, zt3ew(1,1,1,1), imigr, nowe, ml_req1 ) 1485 CALL mpprecv( 2, zt3we(1,1,1,2), imigr, nowe ) 1486 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 1487 END SELECT 1488 ! 1489 ! ! Write lateral conditions 1490 iihom = nlci-nreci 1491 ! 1492 SELECT CASE ( nbondi ) 1493 CASE ( -1 ) 1494 DO jl = 1, jpreci 1495 ptab(iihom+jl,:,:) = ptab(iihom+jl,:,:) + zt3ew(:,jl,:,2) 1496 END DO 1497 CASE ( 0 ) 1498 DO jl = 1, jpreci 1499 ptab(jpreci+jl,:,:) = ptab(jpreci+jl,:,:) + zt3we(:,jl,:,2) 1500 ptab(iihom +jl,:,:) = ptab(iihom +jl,:,:) + zt3ew(:,jl,:,2) 1501 END DO 1502 CASE ( 1 ) 1503 DO jl = 1, jpreci 1504 ptab(jpreci+jl,:,:) = ptab(jpreci+jl,:,:) + zt3we(:,jl,:,2) 1505 END DO 1506 END SELECT 1507 1508 1509 ! 3. North and south directions 1510 ! ----------------------------- 1511 ! always closed : we play only with the neigbours 1512 ! 1513 IF( nbondj /= 2 ) THEN ! Read lateral conditions 1514 ijhom = nlcj-jprecj 1515 DO jl = 1, jprecj 1516 zt3sn(:,jl,:,1) = ptab(:,ijhom+jl,:) ; ptab(:,ijhom+jl,:) = 0.0_wp 1517 zt3ns(:,jl,:,1) = ptab(:,jl ,:) ; ptab(:,jl ,:) = 0.0_wp 1518 END DO 1519 ENDIF 1520 ! 1521 ! ! Migrations 1522 imigr = jprecj * jpi * jpk 1523 ! 1524 SELECT CASE ( nbondj ) 1525 CASE ( -1 ) 1526 CALL mppsend( 4, zt3sn(1,1,1,1), imigr, nono, ml_req1 ) 1527 CALL mpprecv( 3, zt3ns(1,1,1,2), imigr, nono ) 1528 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 1529 CASE ( 0 ) 1530 CALL mppsend( 3, zt3ns(1,1,1,1), imigr, noso, ml_req1 ) 1531 CALL mppsend( 4, zt3sn(1,1,1,1), imigr, nono, ml_req2 ) 1532 CALL mpprecv( 3, zt3ns(1,1,1,2), imigr, nono ) 1533 CALL mpprecv( 4, zt3sn(1,1,1,2), imigr, noso ) 1534 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 1535 IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err) 1536 CASE ( 1 ) 1537 CALL mppsend( 3, zt3ns(1,1,1,1), imigr, noso, ml_req1 ) 1538 CALL mpprecv( 4, zt3sn(1,1,1,2), imigr, noso ) 1539 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 1540 END SELECT 1541 ! 1542 ! ! Write lateral conditions 1543 ijhom = nlcj-nrecj 1544 ! 1545 SELECT CASE ( nbondj ) 1546 CASE ( -1 ) 1547 DO jl = 1, jprecj 1548 ptab(:,ijhom+jl,:) = ptab(:,ijhom+jl,:) + zt3ns(:,jl,:,2) 1549 END DO 1550 CASE ( 0 ) 1551 DO jl = 1, jprecj 1552 ptab(:,jprecj+jl,:) = ptab(:,jprecj+jl,:) + zt3sn(:,jl,:,2) 1553 ptab(:,ijhom +jl,:) = ptab(:,ijhom +jl,:) + zt3ns(:,jl,:,2) 1554 END DO 1555 CASE ( 1 ) 1556 DO jl = 1, jprecj 1557 ptab(:,jprecj+jl,:) = ptab(:,jprecj+jl,:) + zt3sn(:,jl ,:,2) 1558 END DO 1559 END SELECT 1560 1561 1562 ! 4. north fold treatment 1563 ! ----------------------- 1564 ! 1565 IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN 1566 ! 1567 SELECT CASE ( jpni ) 1568 CASE ( 1 ) ; CALL lbc_nfd ( ptab, cd_type, psgn ) ! only 1 northern proc, no mpp 1569 CASE DEFAULT ; CALL mpp_lbc_north( ptab, cd_type, psgn ) ! for all northern procs. 1570 END SELECT 1571 ! 1572 ENDIF 1573 ! 1574 DEALLOCATE( zt3ns, zt3sn, zt3ew, zt3we ) 1575 ! 1576 END SUBROUTINE mpp_lnk_sum_3d 1577 1578 SUBROUTINE mpp_lnk_sum_2d( pt2d, cd_type, psgn, cd_mpp, pval ) 1579 !!---------------------------------------------------------------------- 1580 !! *** routine mpp_lnk_sum_2d *** 1581 !! 1582 !! ** Purpose : Message passing manadgement for 2d array (sum the overlap region) 1583 !! 1584 !! ** Method : Use mppsend and mpprecv function for passing mask 1585 !! between processors following neighboring subdomains. 1586 !! domain parameters 1587 !! nlci : first dimension of the local subdomain 1588 !! nlcj : second dimension of the local subdomain 1589 !! nbondi : mark for "east-west local boundary" 1590 !! nbondj : mark for "north-south local boundary" 1591 !! noea : number for local neighboring processors 1592 !! nowe : number for local neighboring processors 1593 !! noso : number for local neighboring processors 1594 !! nono : number for local neighboring processors 1595 !! 1596 !!---------------------------------------------------------------------- 1597 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pt2d ! 2D array on which the boundary condition is applied 1598 CHARACTER(len=1) , INTENT(in ) :: cd_type ! define the nature of ptab array grid-points 1599 ! ! = T , U , V , F , W and I points 1600 REAL(wp) , INTENT(in ) :: psgn ! =-1 the sign change across the north fold boundary 1601 ! ! = 1. , the sign is kept 1602 CHARACTER(len=3), OPTIONAL , INTENT(in ) :: cd_mpp ! fill the overlap area only 1603 REAL(wp) , OPTIONAL , INTENT(in ) :: pval ! background value (used at closed boundaries) 1604 !! 1605 INTEGER :: ji, jj, jl ! dummy loop indices 1606 INTEGER :: imigr, iihom, ijhom ! temporary integers 1607 INTEGER :: ml_req1, ml_req2, ml_err ! for key_mpi_isend 1608 REAL(wp) :: zland 1609 INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ml_stat ! for key_mpi_isend 1610 ! 1611 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zt2ns, zt2sn ! 2d for north-south & south-north 1612 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zt2ew, zt2we ! 2d for east-west & west-east 1613 1614 !!---------------------------------------------------------------------- 1615 1616 ALLOCATE( zt2ns(jpi,jprecj,2), zt2sn(jpi,jprecj,2), & 1617 & zt2ew(jpj,jpreci,2), zt2we(jpj,jpreci,2) ) 1618 1619 ! 1620 IF( PRESENT( pval ) ) THEN ; zland = pval ! set land value 1621 ELSE ; zland = 0.e0 ! zero by default 1622 ENDIF 1623 1624 ! 1. standard boundary treatment 1625 ! ------------------------------ 1626 ! 2. East and west directions exchange 1627 ! ------------------------------------ 1628 ! we play with the neigbours AND the row number because of the periodicity 1629 ! 1630 SELECT CASE ( nbondi ) ! Read lateral conditions 1631 CASE ( -1, 0, 1 ) ! all exept 2 (i.e. close case) 1632 iihom = nlci - jpreci 1633 DO jl = 1, jpreci 1634 zt2ew(:,jl,1) = pt2d(jl ,:) ; pt2d(jl ,:) = 0.0_wp 1635 zt2we(:,jl,1) = pt2d(iihom +jl,:) ; pt2d(iihom +jl,:) = 0.0_wp 1636 END DO 1637 END SELECT 1638 ! 1639 ! ! Migrations 1640 imigr = jpreci * jpj 1641 ! 1642 SELECT CASE ( nbondi ) 1643 CASE ( -1 ) 1644 CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req1 ) 1645 CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea ) 1646 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1647 CASE ( 0 ) 1648 CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 ) 1649 CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req2 ) 1650 CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea ) 1651 CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe ) 1652 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1653 IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err) 1654 CASE ( 1 ) 1655 CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 ) 1656 CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe ) 1657 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1658 END SELECT 1659 ! 1660 ! ! Write lateral conditions 1661 iihom = nlci-nreci 1662 ! 1663 SELECT CASE ( nbondi ) 1664 CASE ( -1 ) 1665 DO jl = 1, jpreci 1666 pt2d(iihom+jl,:) = pt2d(iihom+jl,:) + zt2ew(:,jl,2) 1667 END DO 1668 CASE ( 0 ) 1669 DO jl = 1, jpreci 1670 pt2d(jpreci+jl,:) = pt2d(jpreci+jl,:) + zt2we(:,jl,2) 1671 pt2d(iihom +jl,:) = pt2d(iihom +jl,:) + zt2ew(:,jl,2) 1672 END DO 1673 CASE ( 1 ) 1674 DO jl = 1, jpreci 1675 pt2d(jpreci+jl,:) = pt2d(jpreci+jl,:) + zt2we(:,jl,2) 1676 END DO 1677 END SELECT 1678 1679 1680 ! 3. North and south directions 1681 ! ----------------------------- 1682 ! always closed : we play only with the neigbours 1683 ! 1684 IF( nbondj /= 2 ) THEN ! Read lateral conditions 1685 ijhom = nlcj - jprecj 1686 DO jl = 1, jprecj 1687 zt2sn(:,jl,1) = pt2d(:,ijhom +jl) ; pt2d(:,ijhom +jl) = 0.0_wp 1688 zt2ns(:,jl,1) = pt2d(:,jl ) ; pt2d(:,jl ) = 0.0_wp 1689 END DO 1690 ENDIF 1691 ! 1692 ! ! Migrations 1693 imigr = jprecj * jpi 1694 ! 1695 SELECT CASE ( nbondj ) 1696 CASE ( -1 ) 1697 CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req1 ) 1698 CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono ) 1699 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1700 CASE ( 0 ) 1701 CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 ) 1702 CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req2 ) 1703 CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono ) 1704 CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso ) 1705 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1706 IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err) 1707 CASE ( 1 ) 1708 CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 ) 1709 CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso ) 1710 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1711 END SELECT 1712 ! 1713 ! ! Write lateral conditions 1714 ijhom = nlcj-nrecj 1715 ! 1716 SELECT CASE ( nbondj ) 1717 CASE ( -1 ) 1718 DO jl = 1, jprecj 1719 pt2d(:,ijhom+jl) = pt2d(:,ijhom+jl) + zt2ns(:,jl,2) 1720 END DO 1721 CASE ( 0 ) 1722 DO jl = 1, jprecj 1723 pt2d(:,jprecj+jl) = pt2d(:,jprecj+jl) + zt2sn(:,jl,2) 1724 pt2d(:,ijhom +jl) = pt2d(:,ijhom +jl) + zt2ns(:,jl,2) 1725 END DO 1726 CASE ( 1 ) 1727 DO jl = 1, jprecj 1728 pt2d(:,jprecj+jl) = pt2d(:,jprecj+jl) + zt2sn(:,jl,2) 1729 END DO 1730 END SELECT 1731 1732 1733 ! 4. north fold treatment 1734 ! ----------------------- 1735 ! 1736 IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN 1737 ! 1738 SELECT CASE ( jpni ) 1739 CASE ( 1 ) ; CALL lbc_nfd ( pt2d, cd_type, psgn ) ! only 1 northern proc, no mpp 1740 CASE DEFAULT ; CALL mpp_lbc_north( pt2d, cd_type, psgn ) ! for all northern procs. 1741 END SELECT 1742 ! 1743 ENDIF 1744 ! 1745 DEALLOCATE( zt2ns, zt2sn, zt2ew, zt2we ) 1746 ! 1747 END SUBROUTINE mpp_lnk_sum_2d 1404 1748 1405 1749 SUBROUTINE mppsend( ktyp, pmess, kbytes, kdest, md_req ) -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90
r5130 r6006 136 136 137 137 imask(:,:)=1 138 WHERE ( zdta(:,:) - zdtaisf(:,:) <= 0.) imask = 0138 WHERE ( zdta(:,:) - zdtaisf(:,:) <= rn_isfhmin ) imask = 0 139 139 140 140 ! 1. Dimension arrays for subdomains -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
r5836 r6006 119 119 IF ( lwp ) WRITE(numout,*) ' ln_divisf = ', ln_divisf 120 120 IF ( lwp ) WRITE(numout,*) ' nn_gammablk = ', nn_gammablk 121 IF ( lwp ) WRITE(numout,*) ' rn_gammat0 = ', rn_gammat0 122 IF ( lwp ) WRITE(numout,*) ' rn_gammas0 = ', rn_gammas0 121 123 IF ( lwp ) WRITE(numout,*) ' rn_tfri2 = ', rn_tfri2 122 124 IF (ln_divisf) THEN ! keep it in the namelist ??? used true anyway as for runoff ? (PM) … … 230 232 CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') 231 233 CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T') 232 CALL sbc_isf_tbl(un(:,:,:) ,utbl(:,:),'U')233 CALL sbc_isf_tbl(vn(:,:,:) ,vtbl(:,:),'V')234 CALL sbc_isf_tbl(un(:,:,:) ,utbl(:,:),'U') 235 CALL sbc_isf_tbl(vn(:,:,:) ,vtbl(:,:),'V') 234 236 ! iom print 235 237 CALL iom_put('ttbl',ttbl(:,:)) 236 238 CALL iom_put('stbl',stbl(:,:)) 237 CALL iom_put('utbl',utbl(:,:) )238 CALL iom_put('vtbl',vtbl(:,:) )239 CALL iom_put('utbl',utbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 240 CALL iom_put('vtbl',vtbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 239 241 ! compute fwf and heat flux 240 242 CALL sbc_isf_cav (kt) … … 356 358 IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 357 359 ! 358 359 ! This test is false only in the very first time step of a run (JMM ???- Initialy build to skip 1rst year of run )360 360 DO ji = 1, jpi 361 361 DO jj = 1, jpj … … 366 366 ! 3. -----------the average temperature between 200m and 600m --------------------- 367 367 DO jk = misfkt(ji,jj),misfkb(ji,jj) 368 ! freezing point temperature at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!)369 ! after verif with UNESCO, wrong sign in BG eq. 2370 368 ! Calculate freezing temperature 371 zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04 372 CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, zpress) 369 CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, fsdept(ji,jj,ik)) 373 370 zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik) ! sum temp 374 371 ENDDO … … 418 415 REAL(wp) :: zfwflx, zhtflx, zhtflx_b 419 416 REAL(wp) :: zgammat, zgammas 420 REAL(wp) :: zeps = -1.e-20_wp !== Local constant initialization ==!417 REAL(wp) :: zeps = 1.e-20_wp !== Local constant initialization ==! 421 418 INTEGER :: ji, jj ! dummy loop indices 422 419 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers … … 503 500 zeps1=rcp*rau0*zgammat 504 501 zeps2=lfusisf*rau0*zgammas 505 zeps3=rhoisf*rcpi*kappa/ risfdep(ji,jj)502 zeps3=rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps) 506 503 zeps4=zlamb2+zlamb3*risfdep(ji,jj) 507 504 zeps6=zeps4-zti(ji,jj) 508 505 zeps7=zeps4-tsurf 509 506 zaqe=zlamb1 * (zeps1 + zeps3) 510 zaqer=0.5/ zaqe507 zaqer=0.5/MIN(zaqe,-zeps) 511 508 zbqe=zeps1*zeps6+zeps3*zeps7-zeps2 512 509 zcqe=zeps2*stbl(ji,jj) … … 517 514 zfrz(ji,jj)=zeps4+zlamb1*zsfrz 518 515 519 ! zfwflx is upward water flux 520 zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / zsfrz)516 ! zfwflx is upward water flux (MAX usefull if kappa = 0 517 zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) ) 521 518 ! zhtflx is upward heat flux (out of ocean) 522 519 ! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
r5836 r6006 100 100 ! 101 101 CALL wrk_alloc( jpi,jpj,jpk,jpts, zts_dta ) 102 !103 102 IF( l_trdtra ) THEN !* Save ta and sa trends 104 103 CALL wrk_alloc( jpi,jpj,jpk,jpts, ztrdts ) -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
r5643 r6006 22 22 USE sbcrnf ! River runoff 23 23 USE sbcisf ! Ice shelf 24 USE iscplini ! Ice sheet coupling 24 25 USE traqsr ! solar radiation penetration 25 26 USE trd_oce ! trends: ocean variables … … 117 118 INTEGER :: ji, jj, jk, jn ! dummy loop indices 118 119 INTEGER :: ikt, ikb 119 INTEGER :: nk_isf120 120 REAL(wp) :: zfact, z1_e3t, zdep 121 121 REAL(wp) :: zalpha, zhk … … 219 219 ! 220 220 IF( nn_isf > 0 ) THEN 221 zfact = 0.5 e0221 zfact = 0.5_wp 222 222 DO jj = 2, jpj 223 223 DO ji = fs_2, fs_jpim1 … … 230 230 ! sign - because fwf sign of evapo (rnf sign of precip) 231 231 DO jk = ikt, ikb - 1 232 ! compute tfreez for the temperature correction (we add water at freezing temperature)233 232 ! compute trend 234 233 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & … … 239 238 240 239 ! level partially include in ice shelf boundary layer 241 ! compute tfreez for the temperature correction (we add water at freezing temperature)242 240 ! compute trend 243 241 tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem) & … … 279 277 ENDIF 280 278 281 IF( l_trdtra ) THEN ! send trends for further diagnostics 279 !---------------------------------------- 280 ! Ice Sheet coupling imbalance correction to have conservation 281 !---------------------------------------- 282 ! 283 IF( ln_iscpl .AND. ln_hsb) THEN ! input of heat and salt due to river runoff 284 DO jk = 1,jpk 285 DO jj = 2, jpj 286 DO ji = fs_2, fs_jpim1 287 zdep = 1._wp / fse3t_n(ji,jj,jk) 288 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) & 289 & * zdep 290 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) & 291 & * zdep 292 END DO 293 END DO 294 END DO 295 ENDIF 296 297 IF( l_trdtra ) THEN ! save the horizontal diffusive trends for further diagnostics 282 298 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 283 299 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90
r4161 r6006 24 24 PRIVATE 25 25 26 PUBLIC glob_sum ! used in many places 27 PUBLIC DDPDD ! also used in closea module 26 PUBLIC glob_sum ! used in many places (masked with tmask_i) 27 PUBLIC glob_sum_full ! used in many places (masked with tmask_h, ie omly over the halos) 28 PUBLIC DDPDD ! also used in closea module 28 29 PUBLIC glob_min, glob_max 29 30 #if defined key_nosignedzero … … 34 35 MODULE PROCEDURE glob_sum_1d, glob_sum_2d, glob_sum_3d, & 35 36 & glob_sum_2d_a, glob_sum_3d_a 37 END INTERFACE 38 INTERFACE glob_sum_full 39 MODULE PROCEDURE glob_sum_full_2d, glob_sum_full_3d 36 40 END INTERFACE 37 41 INTERFACE glob_min … … 156 160 ! 157 161 END FUNCTION glob_sum_3d_a 162 163 FUNCTION glob_sum_full_2d( ptab ) 164 !!---------------------------------------------------------------------- 165 !! *** FUNCTION glob_sum_full_2d *** 166 !! 167 !! ** Purpose : perform a sum in calling DDPDD routine (nomask) 168 !!---------------------------------------------------------------------- 169 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptab 170 REAL(wp) :: glob_sum_full_2d ! global sum 171 !! 172 !!----------------------------------------------------------------------- 173 ! 174 glob_sum_full_2d = SUM( ptab(:,:) * tmask_h(:,:) ) 175 IF( lk_mpp ) CALL mpp_sum( glob_sum_full_2d ) 176 ! 177 END FUNCTION glob_sum_full_2d 178 179 FUNCTION glob_sum_full_3d( ptab ) 180 !!---------------------------------------------------------------------- 181 !! *** FUNCTION glob_sum_full_3d *** 182 !! 183 !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine (nomask) 184 !!---------------------------------------------------------------------- 185 REAL(wp), INTENT(in), DIMENSION(:,:,:) :: ptab 186 REAL(wp) :: glob_sum_full_3d ! global sum 187 !! 188 INTEGER :: ji, jj, jk ! dummy loop indices 189 INTEGER :: ijpk ! local variables: size of ptab 190 !!----------------------------------------------------------------------- 191 ! 192 ijpk = SIZE(ptab,3) 193 ! 194 glob_sum_full_3d = 0.e0 195 DO jk = 1, ijpk 196 glob_sum_full_3d = glob_sum_full_3d + SUM( ptab(:,:,jk) * tmask_h(:,:) ) 197 END DO 198 IF( lk_mpp ) CALL mpp_sum( glob_sum_full_3d ) 199 ! 200 END FUNCTION glob_sum_full_3d 201 158 202 159 203 #else … … 314 358 END FUNCTION glob_sum_3d_a 315 359 360 FUNCTION glob_sum_full_2d( ptab ) 361 !!---------------------------------------------------------------------- 362 !! *** FUNCTION glob_sum_full_2d *** 363 !! 364 !! ** Purpose : perform a sum in calling DDPDD routine 365 !!---------------------------------------------------------------------- 366 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptab 367 REAL(wp) :: glob_sum_full_2d ! global sum (nomask) 368 !! 369 COMPLEX(wp):: ctmp 370 REAL(wp) :: ztmp 371 INTEGER :: ji, jj ! dummy loop indices 372 !!----------------------------------------------------------------------- 373 ! 374 ztmp = 0.e0 375 ctmp = CMPLX( 0.e0, 0.e0, wp ) 376 DO jj = 1, jpj 377 DO ji =1, jpi 378 ztmp = ptab(ji,jj) * tmask_h(ji,jj) 379 CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 380 END DO 381 END DO 382 IF( lk_mpp ) CALL mpp_sum( ctmp ) ! sum over the global domain 383 glob_sum_full_2d = REAL(ctmp,wp) 384 ! 385 END FUNCTION glob_sum_full_2d 386 387 FUNCTION glob_sum_full_3d( ptab ) 388 !!---------------------------------------------------------------------- 389 !! *** FUNCTION glob_sum_full_3d *** 390 !! 391 !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine 392 !!---------------------------------------------------------------------- 393 REAL(wp), INTENT(in), DIMENSION(:,:,:) :: ptab 394 REAL(wp) :: glob_sum_full_3d ! global sum (nomask) 395 !! 396 COMPLEX(wp):: ctmp 397 REAL(wp) :: ztmp 398 INTEGER :: ji, jj, jk ! dummy loop indices 399 INTEGER :: ijpk ! local variables: size of ptab 400 !!----------------------------------------------------------------------- 401 ! 402 ijpk = SIZE(ptab,3) 403 ! 404 ztmp = 0.e0 405 ctmp = CMPLX( 0.e0, 0.e0, wp ) 406 DO jk = 1, ijpk 407 DO jj = 1, jpj 408 DO ji =1, jpi 409 ztmp = ptab(ji,jj,jk) * tmask_h(ji,jj) 410 CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 411 END DO 412 END DO 413 END DO 414 IF( lk_mpp ) CALL mpp_sum( ctmp ) ! sum over the global domain 415 glob_sum_full_3d = REAL(ctmp,wp) 416 ! 417 END FUNCTION glob_sum_full_3d 418 419 420 316 421 #endif 317 422
Note: See TracChangeset
for help on using the changeset viewer.
