Changeset 6006
- Timestamp:
- 2015-12-04T17:56:07+01:00 (8 years ago)
- Location:
- branches/2015/dev_MetOffice_merge_2015
- Files:
-
- 27 edited
- 4 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Chapters/Chap_DOM.tex
r5120 r6006 498 498 Contrary to the horizontal grid, the vertical grid is computed in the code and no 499 499 provision is made for reading it from a file. The only input file is the bathymetry 500 (in meters) (\ifile{bathy\_meter}) 500 (in meters) (\ifile{bathy\_meter}). 501 501 \footnote{N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the 502 502 \ifile{bathy\_meter} file, so that the computation of the number of wet ocean point 503 503 in each water column is by-passed}. 504 If \np{ln\_isfcav}~=~true, an extra file input file describing the ice shelf draft 505 (in meters) (\ifile{isf\_draft\_meter}) is needed and all the location where the isf cavity thinnest 506 than \np{rn\_isfhmin} meters are grounded (ie masked). 507 504 508 After reading the bathymetry, the algorithm for vertical grid definition differs 505 509 between the different options: -
branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Chapters/Chap_SBC.tex
r5120 r6006 996 996 at each relevant depth level, added to the horizontal divergence (\textit{hdivn}) in the subroutine \rou{sbc\_isf\_div} 997 997 (called from \mdl{divcur}). 998 999 \section{ Ice sheet coupling} 1000 \label{SBC_iscpl} 1001 %------------------------------------------namsbc_iscpl---------------------------------------------------- 1002 \namdisplay{namsbc_iscpl} 1003 %-------------------------------------------------------------------------------------------------------- 1004 Ice sheet/ocean coupling is done through file exchange at the restart step. NEMO, at each restart step, 1005 read the bathymetry and ice shelf draft variable in a netcdf file. 1006 If \np{ln\_iscpl = ~true}, the isf draft is assume to be different at each restart step 1007 with potentially some new wet/dry cells due to the ice sheet dynamics/thermodynamics. 1008 The wetting and drying scheme applied on the restart is very simple and described below for the 6 different cases: 1009 \begin{description} 1010 \item[Thin a cell down:] 1011 T/S/ssh are unchanged and U/V in the top cell are corrected to keep the barotropic transport (bt) constant ($bt_b=bt_n$). 1012 \item[Enlarge a cell:] 1013 See case "Thin a cell down" 1014 \item[Dry a cell:] 1015 mask, T/S, U/V and ssh are set to 0. Furthermore, U/V into the water column are modified to satisfy ($bt_b=bt_n$). 1016 \item[Wet a cell:] 1017 mask is set to 1, T/S is extrapolated from neighbours, $ssh_n = ssh_b$ and U/V set to 0. If no neighbours along i,j and k, T/S/U/V and mask are set to 0. 1018 \item[Dry a column:] 1019 mask, T/S, U/V are set to 0 everywhere in the column and ssh set to 0. 1020 \item[Wet a column:] 1021 set mask to 1, T/S is extrapolated from neighbours, ssh is extrapolated from neighbours and U/V set to 0. If no neighbour, T/S/U/V and mask set to 0. 1022 \end{description} 1023 The extrapolation is call \np{nn\_drown} times. It means that if the grounding line retreat by more than \np{nn\_drown} cells between 2 coupling steps, 1024 the code will be unable to fill all the new wet cells properly. The default number is set up for the MISOMIP idealised experiments.\\ 1025 This coupling procedure is able to take into account grounding line and calving front migration. However, it is a non-conservative processe. 1026 This could lead to a trend in heat/salt content and volume. In order to remove the trend and keep the conservation level as close to 0 as possible, 1027 a simple conservation scheme is available with \np{ln\_hsb = ~true}. The heat/salt/vol. gain/loss is diagnose, as well as the location. 1028 Based on what is done on sbcrnf to prescribed a source of heat/salt/vol., the heat/salt/vol. gain/loss is removed/added, 1029 over a period of \np{rn\_fiscpl} time step, into the system. 1030 So after \np{rn\_fiscpl} time step, all the heat/salt/vol. gain/loss due to extrapolation process is canceled.\\ 1031 1032 As the before and now fields are not compatible (modification of the geometry), the restart time step is prescribed to be an euler time step instead of a leap frog and $fields_b = fields_n$. 998 1033 % 999 1034 % ================================================================ -
branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namdom
r4560 r6006 6 6 nn_msh = 0 ! create (=1) a mesh file or not (=0) 7 7 rn_hmin = -3. ! min depth of the ocean (>0) or min number of ocean level (<0) 8 rn_isfhmin = 1.00 ! treshold (m) to discriminate grounding ice to floating ice 8 9 rn_e3zps_min= 20. ! partial step thickness is set larger than the minimum of 9 10 rn_e3zps_rat= 0.1 ! rn_e3zps_min and rn_e3zps_rat*e3t, with 0<rn_e3zps_rat<1 -
branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namrun
r4147 r6006 15 15 cn_ocerst_in = "restart" ! suffix of ocean restart name (input) 16 16 cn_ocerst_out = "restart" ! suffix of ocean restart name (output) 17 ln_iscpl = .false. ! cavity evolution forcing or coupling to ice sheet model (ln_iscpl = T => namsbc_iscpl) 17 18 nn_istate = 0 ! output the initial state (1) or not (0) 18 19 nn_stock = 5475 ! frequency of creation of a restart file (modulo referenced to 1) -
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.