Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/CONFIG/SHARED/field_def.xml (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/CONFIG/SHARED/field_def.xml (revision 5619)
@@ -193,10 +193,16 @@
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/CONFIG/SHARED/namelist_ref (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/CONFIG/SHARED/namelist_ref (revision 5619)
@@ -40,4 +40,5 @@
cn_ocerst_out = "restart" ! suffix of ocean restart name (output)
cn_ocerst_outdir = "." ! directory in which to write output ocean restarts
+ ln_iscpl = .false. ! cavity evolution forcing or coupling to ice sheet model
nn_istate = 0 ! output the initial state (1) or not (0)
ln_rst_list = .false. ! output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F)
@@ -225,4 +226,5 @@
!! namsbc_rnf river runoffs
!! namsbc_isf ice shelf melting/freezing
+!! namsbc_iscpl coupling option between land ice model and ocean
!! namsbc_apr Atmospheric Pressure
!! namsbc_ssr sea surface restoring term (for T and/or S)
@@ -469,4 +471,10 @@
/
!-----------------------------------------------------------------------
+&namsbc_iscpl ! land ice / ocean coupling option
+!-----------------------------------------------------------------------
+ rn_fiscpl = 43800 ! (number of time step) conservation period (maybe should be fix to the coupling frequencey of restart frequency)
+ ln_hfb = .false. ! activate conservation module (conservation exact after a time of rn_fiscpl)
+/
+!-----------------------------------------------------------------------
&namsbc_apr ! Atmospheric pressure used as ocean forcing or in bulk
!-----------------------------------------------------------------------
@@ -802,4 +810,12 @@
rn_smsh = 1. ! Smagorinsky diffusivity: = 0 - use only sheer
rn_aht_m = 2000. ! upper limit or stability criteria for lateral eddy diffusivity (m2/s)
+/
+!-----------------------------------------------------------------------
+&namtra_dmpfile ! tracer: T & S newtonian damping
+!-----------------------------------------------------------------------
+! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask !
+! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename !
+ sn_dmpt = 'resto', -1 ,'Tinit' , .true. , .true. , 'yearly' , '' , '' , ''
+ sn_dmps = 'resto', -1 ,'Sinit' , .true. , .true. , 'yearly' , '' , '' , ''
/
!-----------------------------------------------------------------------
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90 (revision 5619)
@@ -197,6 +197,6 @@
! Elevation
ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*tmask_i(ji,jj)
- ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*umask_i(ji,jj)
- ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*vmask_i(ji,jj)
+ ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj)
+ ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj)
END DO
END DO
@@ -326,6 +326,6 @@
X1= ana_amp(ji,jj,jh,1)
X2=-ana_amp(ji,jj,jh,2)
- out_u(ji,jj, jh) = X1 * umask_i(ji,jj)
- out_u(ji,jj,nb_ana+jh) = X2 * umask_i(ji,jj)
+ out_u(ji,jj, jh) = X1 * ssumask(ji,jj)
+ out_u(ji,jj,nb_ana+jh) = X2 * ssumask(ji,jj)
ENDDO
ENDDO
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90 (revision 5619)
@@ -165,12 +165,12 @@
DO jk = 1, jpkm1
! volume variation (calculated with scale factors)
- zdiff_v2 = zdiff_v2 + glob_sum( surf(:,:) * tmask(:,:,jk) &
- & * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) )
+ zdiff_v2 = zdiff_v2 + glob_sum( surf(:,:) * ( tmask(:,:,jk) &
+ & * fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) )
! heat content variation
- zdiff_hc = zdiff_hc + glob_sum( surf(:,:) * tmask(:,:,jk) &
- & * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - hc_loc_ini(:,:,jk) ) )
+ zdiff_hc = zdiff_hc + glob_sum( surf(:,:) * ( tmask(:,:,jk) &
+ & * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - hc_loc_ini(:,:,jk) ) )
! salt content variation
- zdiff_sc = zdiff_sc + glob_sum( surf(:,:) * tmask(:,:,jk) &
- & * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - sc_loc_ini(:,:,jk) ) )
+ zdiff_sc = zdiff_sc + glob_sum( surf(:,:) * ( tmask(:,:,jk) &
+ & * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - sc_loc_ini(:,:,jk) ) )
ENDDO
@@ -200,5 +200,7 @@
! ENDIF
!!gm end
-
+
+ IF (lwp) PRINT *, 'ISCPL CONS HEAT ', kt, zdiff_hc / zvol_tot, zdiff_sc / zvol_tot
+ IF (lwp) PRINT *, 'ISCPL CONS VOL ', kt, zdiff_v1 * 1.e-9, zdiff_v2 * 1.e-9
IF( lk_vvl ) THEN
@@ -217,5 +219,5 @@
CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp ) ! Heat content variation (1.e20 J)
CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9 ) ! Salt content variation (psu*km3)
- CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9 ) ! volume ssh variation (km3)
+ CALL iom_put( 'bgvolssh' , (zdiff_v1+zdiff_v2) * 1.e-9 ) ! volume ssh variation (km3)
CALL iom_put( 'bgfrcvol' , frc_v * 1.e-9 ) ! vol - surface forcing (km3)
CALL iom_put( 'bgfrctem' , frc_t / zvol_tot ) ! hc - surface forcing (C)
@@ -277,7 +279,7 @@
ssh_ini(:,:) = sshn(:,:) ! initial ssh
DO jk = 1, jpk
- e3t_ini (:,:,jk) = fse3t_n(:,:,jk) ! initial vertical scale factors
- hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk) ! initial heat content
- sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk) ! initial salt content
+ e3t_ini (:,:,jk) = fse3t_n(:,:,jk) * tmask(:,:,jk) ! initial vertical scale factors
+ hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk) * tmask(:,:,jk) ! initial heat content
+ sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk) * tmask(:,:,jk) ! initial salt content
END DO
frc_v = 0._wp ! volume trend due to forcing
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90 (revision 5619)
@@ -43,4 +43,5 @@
INTEGER , PUBLIC :: nn_closea !: =0 suppress closed sea/lake from the ORCA domain or not (=1)
INTEGER , PUBLIC :: nn_euler !: =0 start with forward time step or not (=1)
+ LOGICAL , PUBLIC :: ln_iscpl !: coupling with ice sheet
LOGICAL , PUBLIC :: ln_crs !: Apply grid coarsening to dynamical model output or online passive tracers
@@ -252,13 +253,13 @@
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt !: vertical index of the bottom last T- ocean level
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbku, mbkv !: vertical index of the bottom last U- and W- ocean level
- REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters)
- REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i, umask_i, vmask_i, fmask_i !: interior domain T-point mask
- REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bmask !: land/ocean mask of barotropic stream function
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters)
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i !: interior domain T-point mask
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bmask !: land/ocean mask of barotropic stream function
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfdep !: top first ocean level (ISF)
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mikt, miku, mikv, mikf !: first wet T-, U-, V-, F- ocean level (ISF)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfdep !: Iceshelf draft (ISF)
- REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask !: surface domain T-point mask
-
+
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask, ssumask, ssvmask, ssfmask !: surface mask at T-,U-, V- and F-pts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask !: land/ocean mask at T-, U-, V- and F-pts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask !: land/ocean mask at WT-, WU- and WV-pts
@@ -389,13 +390,14 @@
& hift (jpi,jpj) , hifu (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(8) )
- ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , &
- & tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), &
- & bmask(jpi,jpj) , &
- & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) )
+ ALLOCATE( mbathy(jpi,jpj) , bathy (jpi,jpj) , &
+ & tmask_i(jpi,jpj) , &
+ & ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , &
+ & bmask(jpi,jpj) , &
+ & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj) , STAT=ierr(9) )
! (ISF) Allocation of basic array
ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj), &
& mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) , &
- & mikf(jpi,jpj), ssmask(jpi,jpj), STAT=ierr(10) )
+ & mikf(jpi,jpj), STAT=ierr(10) )
ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk), &
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 (revision 5619)
@@ -111,9 +111,9 @@
END DO
! ! Inverse of the local depth
- hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
- hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
+ hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
+ hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
CALL dom_stp ! time step
- IF( nmsh /= 0 ) CALL dom_wri ! Create a domain file
+ IF( nmsh /= 0 .AND. .NOT. ln_iscpl ) CALL dom_wri ! Create a domain file
IF( .NOT.ln_rstart ) CALL dom_ctl ! Domain control
!
@@ -138,5 +138,5 @@
& nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl, &
& nn_it000, nn_itend , nn_date0 , nn_leapy , nn_istate , nn_stock , &
- & nn_write, ln_dimgnnn, ln_mskland , ln_cfmeta , ln_clobber, nn_chunksz, nn_euler
+ & nn_write, ln_dimgnnn, ln_mskland , ln_cfmeta , ln_clobber, nn_chunksz, nn_euler, ln_iscpl
NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, &
& nn_acc , rn_atfp , rn_rdt , rn_rdtmin , &
@@ -192,4 +192,5 @@
WRITE(numout,*) ' overwrite an existing file ln_clobber = ', ln_clobber
WRITE(numout,*) ' NetCDF chunksize (bytes) nn_chunksz = ', nn_chunksz
+ WRITE(numout,*) ' IS coupling at the restart step ln_iscpl = ', ln_iscpl
ENDIF
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 (revision 5619)
@@ -264,21 +264,21 @@
END DO
END DO
- ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point
+ ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet cell at u point
DO jj = 1, jpjm1
DO ji = 1, fs_jpim1 ! vector loop
- umask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:)))
- vmask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:)))
+ ssumask(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:)))
+ ssvmask(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:)))
END DO
DO ji = 1, jpim1 ! NO vector opt.
- fmask_i(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) &
+ ssfmask(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) &
& * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
END DO
END DO
- CALL lbc_lnk( umask, 'U', 1._wp ) ! Lateral boundary conditions
- CALL lbc_lnk( vmask, 'V', 1._wp )
- CALL lbc_lnk( fmask, 'F', 1._wp )
- CALL lbc_lnk( umask_i, 'U', 1._wp ) ! Lateral boundary conditions
- CALL lbc_lnk( vmask_i, 'V', 1._wp )
- CALL lbc_lnk( fmask_i, 'F', 1._wp )
+ CALL lbc_lnk( umask , 'U', 1._wp ) ! Lateral boundary conditions
+ CALL lbc_lnk( vmask , 'V', 1._wp )
+ CALL lbc_lnk( fmask , 'F', 1._wp )
+ CALL lbc_lnk( ssumask, 'U', 1._wp ) ! Lateral boundary conditions
+ CALL lbc_lnk( ssvmask, 'V', 1._wp )
+ CALL lbc_lnk( ssfmask, 'F', 1._wp )
! 3. Ocean/land mask at wu-, wv- and w points
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90 (revision 5619)
@@ -28,5 +28,5 @@
CONTAINS
- SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid )
+ SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid, kkk )
!!----------------------------------------------------------------------
!! *** ROUTINE dom_ngb ***
@@ -40,7 +40,9 @@
REAL(wp) , INTENT(in ) :: plon, plat ! longitude,latitude of the point
INTEGER , INTENT( out) :: kii, kjj ! i-,j-index of the closes grid point
+ INTEGER , INTENT(in ), OPTIONAL :: kkk
CHARACTER(len=1), INTENT(in ) :: cdgrid ! grid name 'T', 'U', 'V', 'W'
!
INTEGER , DIMENSION(2) :: iloc
+ INTEGER :: jk
REAL(wp) :: zlon, zmini
REAL(wp), POINTER, DIMENSION(:,:) :: zglam, zgphi, zmask, zdist
@@ -52,17 +54,23 @@
!
zmask(:,:) = 0._wp
+ jk = 1
+ IF (PRESENT(kkk)) jk=kkk
SELECT CASE( cdgrid )
- CASE( 'U' ) ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,1)
- CASE( 'V' ) ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,1)
- CASE( 'F' ) ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,1)
- CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,1)
+ CASE( 'U' ) ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,jk)
+ CASE( 'V' ) ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,jk)
+ CASE( 'F' ) ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,jk)
+ CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,jk)
END SELECT
- zlon = MOD( plon + 720., 360. ) ! plon between 0 and 360
- zglam(:,:) = MOD( zglam(:,:) + 720., 360. ) ! glam between 0 and 360
- IF( zlon > 270. ) zlon = zlon - 360. ! zlon between -90 and 270
- IF( zlon < 90. ) WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360. ! glam between -180 and 180
-
- zglam(:,:) = zglam(:,:) - zlon
+ IF (jphgr_msh .NE. 2 .AND. jphgr_msh .NE. 3) THEN
+ zlon = MOD( plon + 720., 360. ) ! plon between 0 and 360
+ zglam(:,:) = MOD( zglam(:,:) + 720., 360. ) ! glam between 0 and 360
+ IF( zlon > 270. ) zlon = zlon - 360. ! zlon between -90 and 270
+ IF( zlon < 90. ) WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360. ! glam between -180 and 180
+ zglam(:,:) = zglam(:,:) - zlon
+ ELSE
+ zglam(:,:) = zglam(:,:) - plon
+ END IF
+!
zgphi(:,:) = zgphi(:,:) - plat
zdist(:,:) = zglam(:,:) * zglam(:,:) + zgphi(:,:) * zgphi(:,:)
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 (revision 5619)
@@ -28,5 +28,5 @@
USE in_out_manager ! I/O manager
USE iom ! I/O manager library
- USE restart ! ocean restart
+ USE restart , ONLY : rst_read_open ! ocean restart
USE lib_mpp ! distributed memory computing library
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
@@ -199,6 +199,6 @@
hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
END DO
- hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) )
- hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) )
+ hur_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1. - ssumask(:,:) )
+ hvr_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1. - ssvmask(:,:) )
! Restoring frequencies for z_tilde coordinate
@@ -545,6 +545,6 @@
END DO
! ! Inverse of the local depth
- hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
- hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
+ hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
+ hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 (revision 5619)
@@ -365,5 +365,5 @@
!! - bathy : meter bathymetry (in meters)
!!----------------------------------------------------------------------
- INTEGER :: ji, jj, jl, jk ! dummy loop indices
+ INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: inum ! temporary logical unit
INTEGER :: ierror ! error flag
@@ -472,4 +472,25 @@
risfdep(:,:)=0.e0
misfdep(:,:)=1
+ !
+ ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code
+ IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN
+ risfdep(:,:)=200.e0
+ misfdep(:,:)=1
+ ij0 = 1 ; ij1 = 40
+ DO jj = mj0(ij0), mj1(ij1)
+ risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp
+ END DO
+ WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp
+ !
+ ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN
+ !
+ risfdep(:,:)=0.e0
+ misfdep(:,:)=1
+ ij0 = 1 ; ij1 = 40
+ DO jj = mj0(ij0), mj1(ij1)
+ risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp
+ END DO
+ WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp
+ END IF
!
DEALLOCATE( idta, zdta )
@@ -529,4 +550,9 @@
WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp
END IF
+ ! set grounded point to 0
+ WHERE (bathy(:,:) .LE. risfdep(:,:)+1e-2 )
+ misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
+ mbathy (:,:) = 0 ; bathy (:,:) = 0._wp
+ END WHERE
!
IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
@@ -952,7 +978,6 @@
REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
- REAL(wp) :: zmax ! Maximum depth
REAL(wp) :: zdiff ! temporary scalar
- REAL(wp) :: zrefdep ! temporary scalar
+ REAL(wp) :: zmax ! temporary scalar
REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt
!!---------------------------------------------------------------------
@@ -993,6 +1018,4 @@
END DO
- IF ( ln_isfcav ) CALL zgr_isf
-
! Scale factors and depth at T- and W-points
DO jk = 1, jpk ! intitialization to the reference z-coordinate
@@ -1002,114 +1025,70 @@
e3w_0 (:,:,jk) = e3w_1d (jk)
END DO
- !
- DO jj = 1, jpj
- DO ji = 1, jpi
- ik = mbathy(ji,jj)
- IF( ik > 0 ) THEN ! ocean point only
- ! max ocean level case
- IF( ik == jpkm1 ) THEN
- zdepwp = bathy(ji,jj)
- ze3tp = bathy(ji,jj) - gdepw_1d(ik)
- ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
- e3t_0(ji,jj,ik ) = ze3tp
- e3t_0(ji,jj,ik+1) = ze3tp
- e3w_0(ji,jj,ik ) = ze3wp
- e3w_0(ji,jj,ik+1) = ze3tp
- gdepw_0(ji,jj,ik+1) = zdepwp
- gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp
- gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
- !
- ELSE ! standard case
- IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
- ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
+
+ ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf
+ IF ( ln_isfcav ) CALL zgr_isf
+
+ ! Scale factors and depth at T- and W-points
+ IF ( .NOT. ln_isfcav ) THEN
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ ik = mbathy(ji,jj)
+ IF( ik > 0 ) THEN ! ocean point only
+ ! max ocean level case
+ IF( ik == jpkm1 ) THEN
+ zdepwp = bathy(ji,jj)
+ ze3tp = bathy(ji,jj) - gdepw_1d(ik)
+ ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
+ e3t_0(ji,jj,ik ) = ze3tp
+ e3t_0(ji,jj,ik+1) = ze3tp
+ e3w_0(ji,jj,ik ) = ze3wp
+ e3w_0(ji,jj,ik+1) = ze3tp
+ gdepw_0(ji,jj,ik+1) = zdepwp
+ gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp
+ gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
+ !
+ ELSE ! standard case
+ IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
+ ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
+ ENDIF
+ !gm Bug? check the gdepw_1d
+ ! ... on ik
+ gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &
+ & * ((gdept_1d( ik ) - gdepw_1d(ik) ) &
+ & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))
+ e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &
+ & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )
+ e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &
+ & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
+ ! ... on ik+1
+ e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
+ e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
+ gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
ENDIF
-!gm Bug? check the gdepw_1d
- ! ... on ik
- gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &
- & * ((gdept_1d( ik ) - gdepw_1d(ik) ) &
- & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))
- e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &
- & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )
- e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &
- & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
- ! ... on ik+1
- e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
- e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
- gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
ENDIF
- ENDIF
- END DO
- END DO
- !
- it = 0
- DO jj = 1, jpj
- DO ji = 1, jpi
- ik = mbathy(ji,jj)
- IF( ik > 0 ) THEN ! ocean point only
- e3tp (ji,jj) = e3t_0(ji,jj,ik)
- e3wp (ji,jj) = e3w_0(ji,jj,ik)
- ! test
- zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )
- IF( zdiff <= 0._wp .AND. lwp ) THEN
- it = it + 1
- WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
- WRITE(numout,*) ' bathy = ', bathy(ji,jj)
- WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
- WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )
+ END DO
+ END DO
+ !
+ it = 0
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ ik = mbathy(ji,jj)
+ IF( ik > 0 ) THEN ! ocean point only
+ e3tp (ji,jj) = e3t_0(ji,jj,ik)
+ e3wp (ji,jj) = e3w_0(ji,jj,ik)
+ ! test
+ zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )
+ IF( zdiff <= 0._wp .AND. lwp ) THEN
+ it = it + 1
+ WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
+ WRITE(numout,*) ' bathy = ', bathy(ji,jj)
+ WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
+ WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )
+ ENDIF
ENDIF
- ENDIF
- END DO
- END DO
- !
- IF ( ln_isfcav ) THEN
- ! (ISF) Definition of e3t, u, v, w for ISF case
- DO jj = 1, jpj
- DO ji = 1, jpi
- ik = misfdep(ji,jj)
- IF( ik > 1 ) THEN ! ice shelf point only
- IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik)
- gdepw_0(ji,jj,ik) = risfdep(ji,jj)
-!gm Bug? check the gdepw_0
- ! ... on ik
- gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) &
- & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) &
- & / ( gdepw_1d(ik+1) - gdepw_1d(ik) )
- e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)
- e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
-
- IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column)
- e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)
- ENDIF
- ! ... on ik / ik-1
- e3w_0 (ji,jj,ik ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))
- e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
-! The next line isn't required and doesn't affect results - included for consistency with bathymetry code
- gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
- ENDIF
- END DO
- END DO
- !
- it = 0
- DO jj = 1, jpj
- DO ji = 1, jpi
- ik = misfdep(ji,jj)
- IF( ik > 1 ) THEN ! ice shelf point only
- e3tp (ji,jj) = e3t_0(ji,jj,ik )
- e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )
- ! test
- zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik )
- IF( zdiff <= 0. .AND. lwp ) THEN
- it = it + 1
- WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
- WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)
- WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
- WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj)
- ENDIF
- ENDIF
- END DO
- END DO
+ END DO
+ END DO
END IF
- ! END (ISF)
-
+ !
! Scale factors and depth at U-, V-, UW and VW-points
DO jk = 1, jpk ! initialisation to z-scale factors
@@ -1119,4 +1098,5 @@
e3vw_0(:,:,jk) = e3w_1d(jk)
END DO
+
DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors
DO jj = 1, jpjm1
@@ -1148,4 +1128,5 @@
CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp )
!
+
DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)
@@ -1255,15 +1236,13 @@
!!----------------------------------------------------------------------
!!
- INTEGER :: ji, jj, jk, jl ! dummy loop indices
- INTEGER :: ik, it ! temporary integers
- INTEGER :: id, jd, nprocd
+ INTEGER :: ji, jj, jl, jk ! dummy loop indices
+ INTEGER :: ik, it ! temporary integers
INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF)
- LOGICAL :: ll_print ! Allow control print for debugging
+ REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t
+ REAL(wp) :: zmax ! Maximum and minimum depth
+ REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar
REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
- REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
- REAL(wp) :: zmax, zmin ! Maximum and minimum depth
+ REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t
REAL(wp) :: zdiff ! temporary scalar
- REAL(wp) :: zrefdep ! temporary scalar
- REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar
REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH)
INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH)
@@ -1280,4 +1259,10 @@
ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level
END WHERE
+
+ ! set grounded point to 0
+ WHERE (bathy(:,:) .LE. risfdep(:,:)+1e-2 )
+ misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
+ mbathy (:,:) = 0 ; bathy (:,:) = 0._wp
+ END WHERE
! Compute misfdep for ocean points (i.e. first wet level)
@@ -1292,4 +1277,11 @@
risfdep(:,:) = 0. ; misfdep(:,:) = 1
END WHERE
+
+ ! remove very shallow ice shelf (less than ~ 10m if 75L)
+ IF ( cp_cfg .NE. "isomip" ) THEN
+ WHERE (risfdep(:,:) < 100 )
+ misfdep = 1; risfdep = 0.0_wp;
+ END WHERE
+ END IF
! 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
@@ -1297,5 +1289,5 @@
! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together
DO jl = 1, 10
- WHERE (bathy(:,:) .EQ. risfdep(:,:) )
+ WHERE (bathy(:,:) .LE. risfdep(:,:)+1e-2 )
misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
mbathy (:,:) = 0 ; bathy (:,:) = 0._wp
@@ -1326,11 +1318,11 @@
! split last cell if possible (only where water column is 2 cell or less)
- DO jk = jpkm1, 1, -1
- zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
- WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)
- mbathy(:,:) = jk
- bathy(:,:) = zmax
- END WHERE
- END DO
+ !DO jk = jpkm1, 1, -1
+ ! zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
+ ! WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)
+ ! mbathy(:,:) = jk
+ ! bathy(:,:) = zmax
+ ! END WHERE
+ !END DO
! split top cell if possible (only where water column is 2 cell or less)
@@ -1350,17 +1342,15 @@
! test bathy
IF (risfdep(ji,jj) .GT. 1) THEN
- zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) &
- & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
- zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &
- & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
+ zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
+ zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN
- IF (zbathydiff .LE. zrisfdepdiff) THEN
- bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
- mbathy(ji,jj)= mbathy(ji,jj) + 1
- ELSE
+! IF (zbathydiff .LE. zrisfdepdiff) THEN
+! bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
+! mbathy(ji,jj)= mbathy(ji,jj) + 1
+! ELSE
risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
misfdep(ji,jj) = misfdep(ji,jj) - 1
- END IF
+! END IF
END IF
END IF
@@ -1374,15 +1364,13 @@
! test bathy
IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
- zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1)&
- & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
- zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &
- & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
- IF (zbathydiff .LE. zrisfdepdiff) THEN
- mbathy(ji,jj) = mbathy(ji,jj) + 1
- bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
- ELSE
+! zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
+! zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
+! IF (zbathydiff .LE. zrisfdepdiff) THEN
+! mbathy(ji,jj) = mbathy(ji,jj) + 1
+! bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
+! ELSE
misfdep(ji,jj)= misfdep(ji,jj) - 1
risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
- END IF
+! END IF
ENDIF
END DO
@@ -1393,17 +1381,13 @@
DO ji = 1, jpim1
IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
- zbathydiff =ABS(bathy(ji,jj ) - (gdepw_1d(mbathy (ji,jj)+1) &
- & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat )))
- zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) &
- & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))
- IF (zbathydiff .LE. zrisfdepdiff) THEN
- mbathy(ji,jj) = mbathy(ji,jj) + 1
- bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) &
- & + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat )
- ELSE
+! zbathydiff =ABS(bathy(ji,jj ) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat )))
+! zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))
+! IF (zbathydiff .LE. zrisfdepdiff) THEN
+! mbathy(ji,jj) = mbathy(ji,jj) + 1
+! bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat )
+! ELSE
misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1
- risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) &
- & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
- END IF
+ risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
+! END IF
ENDIF
END DO
@@ -1424,15 +1408,13 @@
DO ji = 1, jpim1
IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN
- zbathydiff =ABS( bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) &
- & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))
- zrisfdepdiff=ABS(risfdep(ji,jj ) - (gdepw_1d(misfdep(ji,jj ) ) &
- & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat )))
- IF (zbathydiff .LE. zrisfdepdiff) THEN
- mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1
- bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )
- ELSE
+! zbathydiff =ABS( bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))
+! zrisfdepdiff=ABS(risfdep(ji,jj ) - (gdepw_1d(misfdep(ji,jj ) ) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat )))
+! IF (zbathydiff .LE. zrisfdepdiff) THEN
+! mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1
+! bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )
+! ELSE
misfdep(ji,jj) = misfdep(ji,jj) - 1
risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat )
- END IF
+! END IF
ENDIF
END DO
@@ -1455,15 +1437,13 @@
DO ji = 1, jpim1
IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
- zbathydiff =ABS( bathy(ji ,jj) - (gdepw_1d(mbathy (ji,jj)+1) &
- & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat )))
- zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) &
- & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))
- IF (zbathydiff .LE. zrisfdepdiff) THEN
- mbathy(ji,jj) = mbathy(ji,jj) + 1
- bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
- ELSE
+! zbathydiff =ABS( bathy(ji ,jj) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat )))
+! zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))
+! IF (zbathydiff .LE. zrisfdepdiff) THEN
+! mbathy(ji,jj) = mbathy(ji,jj) + 1
+! bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
+! ELSE
misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
- END IF
+! END IF
ENDIF
ENDDO
@@ -1485,17 +1465,13 @@
DO ji = 1, jpim1
IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
- zbathydiff =ABS( bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) &
- & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))
- zrisfdepdiff=ABS(risfdep(ji ,jj) - (gdepw_1d(misfdep(ji ,jj) ) &
- & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat )))
- IF (zbathydiff .LE. zrisfdepdiff) THEN
- mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1
- bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) &
- & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )
- ELSE
+! zbathydiff =ABS( bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))
+! zrisfdepdiff=ABS(risfdep(ji ,jj) - (gdepw_1d(misfdep(ji ,jj) ) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat )))
+! IF (zbathydiff .LE. zrisfdepdiff) THEN
+! mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1
+! bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )
+! ELSE
misfdep(ji,jj) = misfdep(ji ,jj) - 1
- risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) &
- & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat )
- END IF
+ risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat )
+! END IF
ENDIF
ENDDO
@@ -1729,8 +1705,4 @@
! end check compatibility ice shelf/bathy
! remove very shallow ice shelf (less than ~ 10m if 75L)
- WHERE (misfdep(:,:) <= 5)
- misfdep = 1; risfdep = 0.0_wp;
- END WHERE
-
IF( icompt == 0 ) THEN
IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry'
@@ -1738,4 +1710,109 @@
IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'
ENDIF
+
+ ! compute scale factor and depth at T- and W- points
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ ik = mbathy(ji,jj)
+ IF( ik > 0 ) THEN ! ocean point only
+ ! max ocean level case
+ IF( ik == jpkm1 ) THEN
+ zdepwp = bathy(ji,jj)
+ ze3tp = bathy(ji,jj) - gdepw_1d(ik)
+ ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
+ e3t_0(ji,jj,ik ) = ze3tp
+ e3t_0(ji,jj,ik+1) = ze3tp
+ e3w_0(ji,jj,ik ) = ze3wp
+ e3w_0(ji,jj,ik+1) = ze3tp
+ gdepw_0(ji,jj,ik+1) = zdepwp
+ gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp
+ gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
+ !
+ ELSE ! standard case
+ IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
+ ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
+ ENDIF
+ ! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
+!gm Bug? check the gdepw_1d
+ ! ... on ik
+ gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &
+ & * ((gdept_1d( ik ) - gdepw_1d(ik) ) &
+ & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))
+ e3t_0 (ji,jj,ik ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik )
+ e3w_0 (ji,jj,ik ) = gdept_0(ji,jj,ik ) - gdept_1d(ik-1)
+ ! ... on ik+1
+ e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
+ e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
+ ENDIF
+ ENDIF
+ END DO
+ END DO
+ !
+ it = 0
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ ik = mbathy(ji,jj)
+ IF( ik > 0 ) THEN ! ocean point only
+ e3tp (ji,jj) = e3t_0(ji,jj,ik)
+ e3wp (ji,jj) = e3w_0(ji,jj,ik)
+ ! test
+ zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )
+ IF( zdiff <= 0._wp .AND. lwp ) THEN
+ it = it + 1
+ WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
+ WRITE(numout,*) ' bathy = ', bathy(ji,jj)
+ WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
+ WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )
+ ENDIF
+ ENDIF
+ END DO
+ END DO
+ !
+ ! (ISF) Definition of e3t, u, v, w for ISF case
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ ik = misfdep(ji,jj)
+ IF( ik > 1 ) THEN ! ice shelf point only
+ IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik)
+ gdepw_0(ji,jj,ik) = risfdep(ji,jj)
+!gm Bug? check the gdepw_0
+ ! ... on ik
+ gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) &
+ & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) &
+ & / ( gdepw_1d(ik+1) - gdepw_1d(ik) )
+ e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)
+ e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
+
+ IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column)
+ e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)
+ ENDIF
+ ! ... on ik / ik-1
+ e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))
+ e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
+! The next line isn't required and doesn't affect results - included for consistency with bathymetry code
+ gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
+ ENDIF
+ END DO
+ END DO
+
+ it = 0
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ ik = misfdep(ji,jj)
+ IF( ik > 1 ) THEN ! ice shelf point only
+ e3tp (ji,jj) = e3t_0(ji,jj,ik )
+ e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )
+ ! test
+ zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik )
+ IF( zdiff <= 0. .AND. lwp ) THEN
+ it = it + 1
+ WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
+ WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)
+ WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
+ WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj)
+ ENDIF
+ ENDIF
+ END DO
+ END DO
CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep )
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90 (revision 5619)
@@ -34,4 +34,5 @@
TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tsd ! structure of input SST (file informations, fields read)
+ TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tsddmp ! structure of input SST (file informations, fields read)
!! * Substitutions
@@ -60,6 +61,8 @@
TYPE(FLD_N), DIMENSION( jpts) :: slf_i ! array of namelist informations on the fields to read
TYPE(FLD_N) :: sn_tem, sn_sal
+ TYPE(FLD_N) :: sn_dmpt, sn_dmps
!!
NAMELIST/namtsd/ ln_tsd_init, ln_tsd_tradmp, cn_dir, sn_tem, sn_sal
+ NAMELIST/namtra_dmpfile/ sn_dmpt, sn_dmps
INTEGER :: ios
!!----------------------------------------------------------------------
@@ -78,4 +81,12 @@
902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtsd in configuration namelist', lwp )
IF(lwm) WRITE ( numond, namtsd )
+
+ REWIND( numnam_ref ) ! Namelist namtra_dmp in reference namelist : Temperature and salinity damping term
+ READ ( numnam_ref, namtra_dmpfile, IOSTAT = ios, ERR = 903)
+903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in reference namelist', lwp )
+
+ REWIND( numnam_cfg ) ! Namelist namtra_dmp in configuration namelist : Temperature and salinity damping term
+ READ ( numnam_cfg, namtra_dmpfile, IOSTAT = ios, ERR = 904 )
+904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist', lwp )
IF( PRESENT( ld_tradmp ) ) ln_tsd_tradmp = .TRUE. ! forces the initialization when tradmp is used
@@ -105,4 +116,5 @@
!
ALLOCATE( sf_tsd(jpts), STAT=ierr0 )
+ ALLOCATE( sf_tsddmp(jpts), STAT=ierr0 )
IF( ierr0 > 0 ) THEN
CALL ctl_stop( 'dta_tsd_init: unable to allocate sf_tsd structure' ) ; RETURN
@@ -113,4 +125,9 @@
ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk) , STAT=ierr2 )
IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )
+ ! dmp file
+ ALLOCATE( sf_tsddmp(jp_tem)%fnow(jpi,jpj,jpk) , STAT=ierr0 )
+ IF( sn_dmpt%ln_tint ) ALLOCATE( sf_tsddmp(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 )
+ ALLOCATE( sf_tsddmp(jp_sal)%fnow(jpi,jpj,jpk) , STAT=ierr2 )
+ IF( sn_dmps%ln_tint ) ALLOCATE( sf_tsddmp(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )
!
IF( ierr0 + ierr1 + ierr2 + ierr3 > 0 ) THEN
@@ -120,4 +137,6 @@
slf_i(jp_tem) = sn_tem ; slf_i(jp_sal) = sn_sal
CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd' )
+ slf_i(jp_tem) = sn_dmpt ; slf_i(jp_sal) = sn_dmps
+ CALL fld_fill( sf_tsddmp, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd' )
!
ENDIF
@@ -128,5 +147,5 @@
- SUBROUTINE dta_tsd( kt, ptsd )
+ SUBROUTINE dta_tsd( kt, ptsd, ptsddmp )
!!----------------------------------------------------------------------
!! *** ROUTINE dta_tsd ***
@@ -145,4 +164,5 @@
INTEGER , INTENT(in ) :: kt ! ocean time-step
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: ptsd ! T & S data
+ REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), OPTIONAL, INTENT( out) :: ptsddmp ! T & S data
!
INTEGER :: ji, jj, jk, jl, jkk ! dummy loop indicies
@@ -155,4 +175,9 @@
!
CALL fld_read( kt, 1, sf_tsd ) !== read T & S data at kt time step ==!
+ IF ( PRESENT(ptsddmp) ) THEN
+ CALL fld_read( kt, 1, sf_tsddmp ) !== read T & S data at kt time step ==!
+ ptsddmp(:,:,:,jp_tem) = sf_tsddmp(jp_tem)%fnow(:,:,:) ! NO mask
+ ptsddmp(:,:,:,jp_sal) = sf_tsddmp(jp_sal)%fnow(:,:,:)
+ END IF
!
!
@@ -304,4 +329,10 @@
IF( sf_tsd(jp_sal)%ln_tint ) DEALLOCATE( sf_tsd(jp_sal)%fdta )
DEALLOCATE( sf_tsd ) ! the structure itself
+ IF(lwp) WRITE(numout,*) 'dta_tsd: deallocte T & S arrays as they are only use to initialize the run'
+ DEALLOCATE( sf_tsddmp(jp_tem)%fnow ) ! T arrays in the structure
+ IF( sf_tsddmp(jp_tem)%ln_tint ) DEALLOCATE( sf_tsddmp(jp_tem)%fdta )
+ DEALLOCATE( sf_tsddmp(jp_sal)%fnow ) ! S arrays in the structure
+ IF( sf_tsddmp(jp_sal)%ln_tint ) DEALLOCATE( sf_tsddmp(jp_sal)%fdta )
+ DEALLOCATE( sf_tsddmp ) ! the structure itself
ENDIF
!
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90 (revision 5619)
@@ -47,4 +47,5 @@
USE wrk_nemo ! Memory allocation
USE timing ! Timing
+ USE sbc_iscpl
IMPLICIT NONE
@@ -91,4 +92,5 @@
! ! -------------------
CALL rst_read ! Read the restart file
+ IF (ln_iscpl) CALL rst_iscpl ! extraloate restart to wet and dry
CALL day_init ! model calendar (using both namelist and restart infos)
ELSE
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90 (revision 5619)
@@ -73,5 +73,5 @@
REAL(wp), PUBLIC :: xlsn !: = lfus*rhosn (volumetric latent heat fusion of snow) [J/m3]
#else
- REAL(wp), PUBLIC :: rhoic = 900._wp !: volumic mass of sea ice [kg/m3]
+ REAL(wp), PUBLIC :: rhoic = 917._wp !: volumic mass of sea ice [kg/m3]
REAL(wp), PUBLIC :: rcdic = 2.034396_wp !: conductivity of the ice [W/m/K]
REAL(wp), PUBLIC :: rcpic = 1.8837e+6_wp !: volumetric specific heat for ice [J/m3/K]
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90 (revision 5619)
@@ -29,4 +29,5 @@
USE sbcrnf ! river runoff
USE sbcisf ! ice shelf
+ USE sbc_iscpl ! ice shelf
USE cla ! cross land advection (cla_div routine)
USE in_out_manager ! I/O manager
@@ -328,5 +329,6 @@
IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field)
- IF( ln_divisf .AND. (nn_isf .GT. 0) ) CALL sbc_isf_div( hdivn ) ! ice shelf (update hdivn field)
+ IF( ln_divisf .AND. (nn_isf .GT. 0) ) CALL sbc_isf_div ( hdivn ) ! ice shelf (update hdivn field)
+ IF( ln_iscpl .AND. ln_hfb ) CALL sbc_iscpl_div( hdivn ) ! ice shelf (update hdivn field)
IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field)
!
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 (revision 5619)
@@ -353,6 +353,6 @@
hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
END DO
- hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) )
- hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) )
+ hur_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
+ hvr_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
ENDIF
!
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90 (revision 5619)
@@ -24,4 +24,5 @@
USE trdmxl_oce ! ocean active mixed layer tracers trends variables
USE divcur ! hor. divergence and curl (div & cur routines)
+ USE wrk_nemo
IMPLICIT NONE
@@ -145,4 +146,17 @@
CALL iom_rstput( kt, nitrst, numrow, 'rhd' , rhd )
#endif
+
+
+
+ IF ( ln_iscpl ) THEN
+ CALL iom_rstput( kt, nitrst, numrow, 'tmask' , tmask ) ! need to extrapolate T/S
+ CALL iom_rstput( kt, nitrst, numrow, 'umask' , umask ) ! need to correct barotropic velocity
+ CALL iom_rstput( kt, nitrst, numrow, 'vmask' , vmask ) ! need to correct barotropic velocity
+ CALL iom_rstput( kt, nitrst, numrow, 'smask' , ssmask ) ! need to correct barotropic velocity
+ CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) ) ! need to compute temperature correction
+ CALL iom_rstput( kt, nitrst, numrow, 'fse3u_n', fse3u_n(:,:,:) ) ! need to compute volume correction ????
+ CALL iom_rstput( kt, nitrst, numrow, 'fse3v_n', fse3v_n(:,:,:) ) ! need to compute volume correction ????
+ CALL iom_rstput( kt, nitrst, numrow, 'fsdepw_n', fsdepw_n(:,:,:) ) ! need to compute volume correction ????
+ END IF
IF( kt == nitrst ) THEN
CALL iom_close( numrow ) ! close the restart file (only at last time step)
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90 (revision 5619)
@@ -31,4 +31,8 @@
END INTERFACE
+ INTERFACE lbc_sum
+ MODULE PROCEDURE mpp_sum_3d, mpp_sum_2d
+ END INTERFACE
+
INTERFACE lbc_bdy_lnk
MODULE PROCEDURE mpp_lnk_bdy_2d, mpp_lnk_bdy_3d
@@ -45,4 +49,5 @@
PUBLIC lbc_lnk ! ocean lateral boundary conditions
PUBLIC lbc_lnk_multi ! modified ocean lateral boundary conditions
+ PUBLIC lbc_sum
PUBLIC lbc_lnk_e
PUBLIC lbc_bdy_lnk ! ocean lateral BDY boundary conditions
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90 (revision 5619)
@@ -72,4 +72,5 @@
PUBLIC mpp_lnk_3d, mpp_lnk_3d_gather, mpp_lnk_2d, mpp_lnk_2d_e
PUBLIC mpp_lnk_2d_9
+ PUBLIC mpp_sum_3d, mpp_sum_2d
PUBLIC mppscatter, mppgather
PUBLIC mpp_ini_ice, mpp_ini_znl
@@ -1394,4 +1395,412 @@
END SUBROUTINE mpp_lnk_2d_e
+ SUBROUTINE mpp_sum_3d( ptab, cd_type, psgn, cd_mpp, pval )
+ !!----------------------------------------------------------------------
+ !! *** routine mpp_sum_3d ***
+ !!
+ !! ** Purpose : Message passing manadgement
+ !!
+ !! ** Method : Use mppsend and mpprecv function for passing mask
+ !! between processors following neighboring subdomains.
+ !! domain parameters
+ !! nlci : first dimension of the local subdomain
+ !! nlcj : second dimension of the local subdomain
+ !! nbondi : mark for "east-west local boundary"
+ !! nbondj : mark for "north-south local boundary"
+ !! noea : number for local neighboring processors
+ !! nowe : number for local neighboring processors
+ !! noso : number for local neighboring processors
+ !! nono : number for local neighboring processors
+ !!
+ !! ** Action : ptab with update value at its periphery
+ !!
+ !!----------------------------------------------------------------------
+ REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: ptab ! 3D array on which the boundary condition is applied
+ CHARACTER(len=1) , INTENT(in ) :: cd_type ! define the nature of ptab array grid-points
+ ! ! = T , U , V , F , W points
+ REAL(wp) , INTENT(in ) :: psgn ! =-1 the sign change across the north fold boundary
+ ! ! = 1. , the sign is kept
+ CHARACTER(len=3), OPTIONAL , INTENT(in ) :: cd_mpp ! fill the overlap area only
+ REAL(wp) , OPTIONAL , INTENT(in ) :: pval ! background value (used at closed boundaries)
+ !!
+ INTEGER :: ji, jj, jk, jl ! dummy loop indices
+ INTEGER :: imigr, iihom, ijhom ! temporary integers
+ INTEGER :: ml_req1, ml_req2, ml_err ! for key_mpi_isend
+ REAL(wp) :: zland
+ INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ml_stat ! for key_mpi_isend
+ !
+ REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zt3ns, zt3sn ! 3d for north-south & south-north
+ REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zt3ew, zt3we ! 3d for east-west & west-east
+
+ !!----------------------------------------------------------------------
+
+ ALLOCATE( zt3ns(jpi,jprecj,jpk,2), zt3sn(jpi,jprecj,jpk,2), &
+ & zt3ew(jpj,jpreci,jpk,2), zt3we(jpj,jpreci,jpk,2) )
+
+ !
+ IF( PRESENT( pval ) ) THEN ; zland = pval ! set land value
+ ELSE ; zland = 0.e0 ! zero by default
+ ENDIF
+
+ ! 1. standard boundary treatment
+ ! ------------------------------
+ IF( PRESENT( cd_mpp ) ) THEN ! only fill added line/raw with existing values
+ !
+ ! WARNING ptab is defined only between nld and nle
+! DO jk = 1, jpk
+! DO jj = nlcj+1, jpj ! added line(s) (inner only)
+! ptab(nldi :nlei , jj ,jk) = ptab(nldi:nlei, nlej,jk)
+! ptab(1 :nldi-1, jj ,jk) = ptab(nldi , nlej,jk)
+! ptab(nlei+1:nlci , jj ,jk) = ptab( nlei, nlej,jk)
+! END DO
+! DO ji = nlci+1, jpi ! added column(s) (full)
+! ptab(ji ,nldj :nlej ,jk) = ptab( nlei,nldj:nlej,jk)
+! ptab(ji ,1 :nldj-1,jk) = ptab( nlei,nldj ,jk)
+! ptab(ji ,nlej+1:jpj ,jk) = ptab( nlei, nlej,jk)
+! END DO
+! END DO
+ !
+ ELSE ! standard close or cyclic treatment
+ !
+ ! ! East-West boundaries
+ ! !* Cyclic east-west
+ IF( nbondi == 2 .AND. (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) ) THEN
+! ptab( 1 ,:,:) = ptab(jpim1,:,:)
+! ptab(jpi,:,:) = ptab( 2 ,:,:)
+ ELSE !* closed
+! IF( .NOT. cd_type == 'F' ) ptab( 1 :jpreci,:,:) = zland ! south except F-point
+! ptab(nlci-jpreci+1:jpi ,:,:) = zland ! north
+ ENDIF
+ ! ! North-South boundaries (always closed)
+! IF( .NOT. cd_type == 'F' ) ptab(:, 1 :jprecj,:) = zland ! south except F-point
+! ptab(:,nlcj-jprecj+1:jpj ,:) = zland ! north
+ !
+ ENDIF
+
+ ! 2. East and west directions exchange
+ ! ------------------------------------
+ ! we play with the neigbours AND the row number because of the periodicity
+ !
+ SELECT CASE ( nbondi ) ! Read Dirichlet lateral conditions
+ CASE ( -1, 0, 1 ) ! all exept 2 (i.e. close case)
+ iihom = nlci-jpreci
+ DO jl = 1, jpreci
+ zt3ew(:,jl,:,1) = ptab(jl ,:,:) ; ptab(jl ,:,:) = 0.0_wp
+ zt3we(:,jl,:,1) = ptab(iihom+jl,:,:) ; ptab(iihom+jl,:,:) = 0.0_wp
+ END DO
+ END SELECT
+ !
+ ! ! Migrations
+ imigr = jpreci * jpj * jpk
+ !
+ SELECT CASE ( nbondi )
+ CASE ( -1 )
+ CALL mppsend( 2, zt3we(1,1,1,1), imigr, noea, ml_req1 )
+ CALL mpprecv( 1, zt3ew(1,1,1,2), imigr, noea )
+ IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
+ CASE ( 0 )
+ CALL mppsend( 1, zt3ew(1,1,1,1), imigr, nowe, ml_req1 )
+ CALL mppsend( 2, zt3we(1,1,1,1), imigr, noea, ml_req2 )
+ CALL mpprecv( 1, zt3ew(1,1,1,2), imigr, noea )
+ CALL mpprecv( 2, zt3we(1,1,1,2), imigr, nowe )
+ IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
+ IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err)
+ CASE ( 1 )
+ CALL mppsend( 1, zt3ew(1,1,1,1), imigr, nowe, ml_req1 )
+ CALL mpprecv( 2, zt3we(1,1,1,2), imigr, nowe )
+ IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
+ END SELECT
+ !
+ ! ! Write Dirichlet lateral conditions
+ iihom = nlci-nreci
+ !
+ SELECT CASE ( nbondi )
+ CASE ( -1 )
+ DO jl = 1, jpreci
+ ptab(iihom+jl,:,:) = ptab(iihom+jl,:,:) + zt3ew(:,jl,:,2)
+ END DO
+ CASE ( 0 )
+ DO jl = 1, jpreci
+ ptab(jpreci+jl,:,:) = ptab(jpreci+jl,:,:) + zt3we(:,jl,:,2)
+ ptab(iihom +jl,:,:) = ptab(iihom +jl,:,:) + zt3ew(:,jl,:,2)
+ END DO
+ CASE ( 1 )
+ DO jl = 1, jpreci
+ ptab(jpreci+jl,:,:) = ptab(jpreci+jl,:,:) + zt3we(:,jl,:,2)
+ END DO
+ END SELECT
+
+
+ ! 3. North and south directions
+ ! -----------------------------
+ ! always closed : we play only with the neigbours
+ !
+ IF( nbondj /= 2 ) THEN ! Read Dirichlet lateral conditions
+ ijhom = nlcj-jprecj
+ DO jl = 1, jprecj
+ zt3sn(:,jl,:,1) = ptab(:,ijhom+jl,:) ; ptab(:,ijhom+jl,:) = 0.0_wp
+ zt3ns(:,jl,:,1) = ptab(:,jl ,:) ; ptab(:,jl ,:) = 0.0_wp
+ END DO
+ ENDIF
+ !
+ ! ! Migrations
+ imigr = jprecj * jpi * jpk
+ !
+ SELECT CASE ( nbondj )
+ CASE ( -1 )
+ CALL mppsend( 4, zt3sn(1,1,1,1), imigr, nono, ml_req1 )
+ CALL mpprecv( 3, zt3ns(1,1,1,2), imigr, nono )
+ IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
+ CASE ( 0 )
+ CALL mppsend( 3, zt3ns(1,1,1,1), imigr, noso, ml_req1 )
+ CALL mppsend( 4, zt3sn(1,1,1,1), imigr, nono, ml_req2 )
+ CALL mpprecv( 3, zt3ns(1,1,1,2), imigr, nono )
+ CALL mpprecv( 4, zt3sn(1,1,1,2), imigr, noso )
+ IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
+ IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err)
+ CASE ( 1 )
+ CALL mppsend( 3, zt3ns(1,1,1,1), imigr, noso, ml_req1 )
+ CALL mpprecv( 4, zt3sn(1,1,1,2), imigr, noso )
+ IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
+ END SELECT
+ !
+ ! ! Write Dirichlet lateral conditions
+ ijhom = nlcj-nrecj
+ !
+ SELECT CASE ( nbondj )
+ CASE ( -1 )
+ DO jl = 1, jprecj
+ ptab(:,ijhom+jl,:) = ptab(:,ijhom+jl,:) + zt3ns(:,jl,:,2)
+ END DO
+ CASE ( 0 )
+ DO jl = 1, jprecj
+ ptab(:,jprecj+jl,:) = ptab(:,jprecj+jl,:) + zt3sn(:,jl,:,2)
+ ptab(:,ijhom +jl,:) = ptab(:,ijhom +jl,:) + zt3ns(:,jl,:,2)
+ END DO
+ CASE ( 1 )
+ DO jl = 1, jprecj
+ ptab(:,jprecj+jl,:) = ptab(:,jprecj+jl,:) + zt3sn(:,jl ,:,2)
+ END DO
+ END SELECT
+
+
+ ! 4. north fold treatment
+ ! -----------------------
+ !
+ IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN
+ !
+ SELECT CASE ( jpni )
+ CASE ( 1 ) ; CALL lbc_nfd ( ptab, cd_type, psgn ) ! only 1 northern proc, no mpp
+ CASE DEFAULT ; CALL mpp_lbc_north( ptab, cd_type, psgn ) ! for all northern procs.
+ END SELECT
+ !
+ ENDIF
+ !
+ DEALLOCATE( zt3ns, zt3sn, zt3ew, zt3we )
+ !
+ END SUBROUTINE mpp_sum_3d
+
+ SUBROUTINE mpp_sum_2d( pt2d, cd_type, psgn, cd_mpp, pval )
+ !!----------------------------------------------------------------------
+ !! *** routine mpp_sum_2d ***
+ !!
+ !! ** Purpose : Message passing manadgement for 2d array
+ !!
+ !! ** Method : Use mppsend and mpprecv function for passing mask
+ !! between processors following neighboring subdomains.
+ !! domain parameters
+ !! nlci : first dimension of the local subdomain
+ !! nlcj : second dimension of the local subdomain
+ !! nbondi : mark for "east-west local boundary"
+ !! nbondj : mark for "north-south local boundary"
+ !! noea : number for local neighboring processors
+ !! nowe : number for local neighboring processors
+ !! noso : number for local neighboring processors
+ !! nono : number for local neighboring processors
+ !!
+ !!----------------------------------------------------------------------
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pt2d ! 2D array on which the boundary condition is applied
+ CHARACTER(len=1) , INTENT(in ) :: cd_type ! define the nature of ptab array grid-points
+ ! ! = T , U , V , F , W and I points
+ REAL(wp) , INTENT(in ) :: psgn ! =-1 the sign change across the north fold boundary
+ ! ! = 1. , the sign is kept
+ CHARACTER(len=3), OPTIONAL , INTENT(in ) :: cd_mpp ! fill the overlap area only
+ REAL(wp) , OPTIONAL , INTENT(in ) :: pval ! background value (used at closed boundaries)
+ !!
+ INTEGER :: ji, jj, jl ! dummy loop indices
+ INTEGER :: imigr, iihom, ijhom ! temporary integers
+ INTEGER :: ml_req1, ml_req2, ml_err ! for key_mpi_isend
+ REAL(wp) :: zland
+ INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ml_stat ! for key_mpi_isend
+ !
+ REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zt2ns, zt2sn ! 2d for north-south & south-north
+ REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zt2ew, zt2we ! 2d for east-west & west-east
+
+ !!----------------------------------------------------------------------
+
+ ALLOCATE( zt2ns(jpi,jprecj,2), zt2sn(jpi,jprecj,2), &
+ & zt2ew(jpj,jpreci,2), zt2we(jpj,jpreci,2) )
+
+ !
+ IF( PRESENT( pval ) ) THEN ; zland = pval ! set land value
+ ELSE ; zland = 0.e0 ! zero by default
+ ENDIF
+
+ ! 1. standard boundary treatment
+ ! ------------------------------
+ !
+! IF( PRESENT( cd_mpp ) ) THEN ! only fill added line/raw with existing values
+! !
+! ! WARNING pt2d is defined only between nld and nle
+! DO jj = nlcj+1, jpj ! added line(s) (inner only)
+! pt2d(nldi :nlei , jj ) = pt2d(nldi:nlei, nlej)
+! pt2d(1 :nldi-1, jj ) = pt2d(nldi , nlej)
+! pt2d(nlei+1:nlci , jj ) = pt2d( nlei, nlej)
+! END DO
+! DO ji = nlci+1, jpi ! added column(s) (full)
+! pt2d(ji ,nldj :nlej ) = pt2d( nlei,nldj:nlej)
+! pt2d(ji ,1 :nldj-1) = pt2d( nlei,nldj )
+! pt2d(ji ,nlej+1:jpj ) = pt2d( nlei, nlej)
+! END DO
+! !
+! ELSE ! standard close or cyclic treatment
+! !
+! ! ! East-West boundaries
+! IF( nbondi == 2 .AND. & ! Cyclic east-west
+! & (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) ) THEN
+! pt2d( 1 ,:) = pt2d(jpim1,:) ! west
+! pt2d(jpi,:) = pt2d( 2 ,:) ! east
+! ELSE ! closed
+! IF( .NOT. cd_type == 'F' ) pt2d( 1 :jpreci,:) = zland ! south except F-point
+! pt2d(nlci-jpreci+1:jpi ,:) = zland ! north
+! ENDIF
+! ! ! North-South boundaries (always closed)
+! IF( .NOT. cd_type == 'F' ) pt2d(:, 1 :jprecj) = zland !south except F-point
+! pt2d(:,nlcj-jprecj+1:jpj ) = zland ! north
+! !
+! ENDIF
+
+ ! 2. East and west directions exchange
+ ! ------------------------------------
+ ! we play with the neigbours AND the row number because of the periodicity
+ !
+ SELECT CASE ( nbondi ) ! Read Dirichlet lateral conditions
+ CASE ( -1, 0, 1 ) ! all exept 2 (i.e. close case)
+ iihom = nlci - jpreci
+ DO jl = 1, jpreci
+ zt2ew(:,jl,1) = pt2d(jl ,:) ; pt2d(jl ,:) = 0.0_wp
+ zt2we(:,jl,1) = pt2d(iihom +jl,:) ; pt2d(iihom +jl,:) = 0.0_wp
+ END DO
+ END SELECT
+ !
+ ! ! Migrations
+ imigr = jpreci * jpj
+ !
+ SELECT CASE ( nbondi )
+ CASE ( -1 )
+ CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req1 )
+ CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea )
+ IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+ CASE ( 0 )
+ CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 )
+ CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req2 )
+ CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea )
+ CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe )
+ IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+ IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
+ CASE ( 1 )
+ CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 )
+ CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe )
+ IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+ END SELECT
+ !
+ ! ! Write Dirichlet lateral conditions
+ iihom = nlci-nreci
+ !
+ SELECT CASE ( nbondi )
+ CASE ( -1 )
+ DO jl = 1, jpreci
+ pt2d(iihom+jl,:) = pt2d(iihom+jl,:) + zt2ew(:,jl,2)
+ END DO
+ CASE ( 0 )
+ DO jl = 1, jpreci
+ pt2d(jpreci+jl,:) = pt2d(jpreci+jl,:) + zt2we(:,jl,2)
+ pt2d(iihom +jl,:) = pt2d(iihom +jl,:) + zt2ew(:,jl,2)
+ END DO
+ CASE ( 1 )
+ DO jl = 1, jpreci
+ pt2d(jpreci+jl,:) = pt2d(jpreci+jl,:) + zt2we(:,jl,2)
+ END DO
+ END SELECT
+
+
+ ! 3. North and south directions
+ ! -----------------------------
+ ! always closed : we play only with the neigbours
+ !
+ IF( nbondj /= 2 ) THEN ! Read Dirichlet lateral conditions
+ ijhom = nlcj - jprecj
+ DO jl = 1, jprecj
+ zt2sn(:,jl,1) = pt2d(:,ijhom +jl) ; pt2d(:,ijhom +jl) = 0.0_wp
+ zt2ns(:,jl,1) = pt2d(:,jl ) ; pt2d(:,jl ) = 0.0_wp
+ END DO
+ ENDIF
+ !
+ ! ! Migrations
+ imigr = jprecj * jpi
+ !
+ SELECT CASE ( nbondj )
+ CASE ( -1 )
+ CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req1 )
+ CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono )
+ IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+ CASE ( 0 )
+ CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 )
+ CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req2 )
+ CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono )
+ CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso )
+ IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+ IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
+ CASE ( 1 )
+ CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 )
+ CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso )
+ IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+ END SELECT
+ !
+ ! ! Write Dirichlet lateral conditions
+ ijhom = nlcj-nrecj
+ !
+ SELECT CASE ( nbondj )
+ CASE ( -1 )
+ DO jl = 1, jprecj
+ pt2d(:,ijhom+jl) = pt2d(:,ijhom+jl) + zt2ns(:,jl,2)
+ END DO
+ CASE ( 0 )
+ DO jl = 1, jprecj
+ pt2d(:,jprecj+jl) = pt2d(:,jprecj+jl) + zt2sn(:,jl,2)
+ pt2d(:,ijhom +jl) = pt2d(:,ijhom +jl) + zt2ns(:,jl,2)
+ END DO
+ CASE ( 1 )
+ DO jl = 1, jprecj
+ pt2d(:,jprecj+jl) = pt2d(:,jprecj+jl) + zt2sn(:,jl,2)
+ END DO
+ END SELECT
+
+
+ ! 4. north fold treatment
+ ! -----------------------
+ !
+ IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN
+ !
+ SELECT CASE ( jpni )
+ CASE ( 1 ) ; CALL lbc_nfd ( pt2d, cd_type, psgn ) ! only 1 northern proc, no mpp
+ CASE DEFAULT ; CALL mpp_lbc_north( pt2d, cd_type, psgn ) ! for all northern procs.
+ END SELECT
+ !
+ ENDIF
+ !
+ DEALLOCATE( zt2ns, zt2sn, zt2ew, zt2we )
+ !
+ END SUBROUTINE mpp_sum_2d
SUBROUTINE mppsend( ktyp, pmess, kbytes, kdest, md_req )
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90 (revision 5619)
@@ -136,5 +136,5 @@
imask(:,:)=1
- WHERE ( zdta(:,:) - zdtaisf(:,:) <= 0. ) imask = 0
+ WHERE ( zdta(:,:) - zdtaisf(:,:) <= 1e-2 ) imask = 0
! 1. Dimension arrays for subdomains
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90 (revision 5619)
@@ -111,4 +111,5 @@
zcoef = z_fwf * rcp
emp(:,:) = emp(:,:) - z_fwf * tmask(:,:,1)
+ sfx(:,:) = sfx(:,:) + z_fwf * sss_m * tmask(:,:,1)
qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
ENDIF
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbciscpl.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbciscpl.F90 (revision 5619)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbciscpl.F90 (revision 5619)
@@ -0,0 +1,811 @@
+MODULE sbc_iscpl
+ !!======================================================================
+ !! *** MODULE sbciscpl***
+ !! Ocean forcing: river runoff
+ !!=====================================================================
+ !! History : NEMO ! 2015-01 P. Mathiot: original
+ !!----------------------------------------------------------------------
+
+ !!----------------------------------------------------------------------
+ !! sbc_iscpl_alloc: variable allocation
+ !! rst_iscpl : restart correction in case of coupling with ice sheet
+ !! sbc_iscpl_div : correction of divergence to keep volume conservation
+ !!----------------------------------------------------------------------
+ USE dom_oce ! ocean space and time domain
+ USE domwri ! ocean space and time domain
+ USE domvvl, ONLY : dom_vvl_interpol
+ USE phycst ! physical constants
+ USE sbc_oce ! surface boundary condition variables
+ USE oce ! global tra/dyn variable
+ USE in_out_manager ! I/O manager
+ USE iom ! I/O module
+ USE lib_mpp ! MPP library
+ USE lib_fortran ! MPP library
+ USE wrk_nemo ! Memory allocation
+ USE lbclnk
+ USE domngb
+ USE iom
+ USE sbc_ice, ONLY : lk_lim3
+
+ IMPLICIT NONE
+ PRIVATE
+
+ PUBLIC sbc_iscpl_div
+ PUBLIC rst_iscpl
+ PUBLIC rst_iscpl_interpol ! routine to wet and dry
+ PUBLIC rst_iscpl_alloc
+ !! !!* namsbc_iscpl namelist *
+ LOGICAL , PUBLIC :: ln_hfb
+ REAL(wp), PUBLIC :: rn_fiscpl
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,: ) :: hdiv_iscpl
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: htsc_iscpl
+ !! * Substitutions
+# include "domzgr_substitute.h90"
+ !!----------------------------------------------------------------------
+ !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+ !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
+ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+ !!----------------------------------------------------------------------
+CONTAINS
+
+ INTEGER FUNCTION rst_iscpl_alloc()
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE sbc_iscpl_alloc ***
+ !!----------------------------------------------------------------------
+ ALLOCATE( htsc_iscpl(jpi,jpj,jpk,jpts) , hdiv_iscpl(jpi,jpj,jpk) , STAT=rst_iscpl_alloc )
+ !
+ IF( lk_mpp ) CALL mpp_sum ( rst_iscpl_alloc )
+ IF( rst_iscpl_alloc > 0 ) CALL ctl_warn('rst_iscpl_alloc: allocation of arrays failed')
+ END FUNCTION rst_iscpl_alloc
+
+ SUBROUTINE rst_iscpl
+ REAL(wp), DIMENSION(:,: ), POINTER :: zsmask_b
+ REAL(wp), DIMENSION(:,:,:), POINTER :: ztmask_b, zumask_b, zvmask_b
+ REAL(wp), DIMENSION(:,:,:), POINTER :: ze3t_b , ze3u_b , ze3v_b
+ REAL(wp), DIMENSION(:,:,:), POINTER :: zdepw_b
+ INTEGER :: jk
+ LOGICAL :: llok
+!!----------------------------------------------------------------------
+ INTEGER :: inum0
+ CHARACTER(20) :: cfile
+
+
+ WRITE(numout,*) 'rst_read ln_iscpl : ', ln_iscpl
+ IF ( ln_iscpl ) THEN
+ CALL wrk_alloc(jpi,jpj,jpk, ztmask_b, zumask_b, zvmask_b) ! need this new variable for evoving isf cavity
+ CALL wrk_alloc(jpi,jpj,jpk, ze3t_b , ze3u_b , ze3v_b ) ! need of this variable afterwards (only used for interpolation)
+ CALL wrk_alloc(jpi,jpj,jpk, zdepw_b )
+ CALL wrk_alloc(jpi,jpj, zsmask_b )
+
+ IF( rst_iscpl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' )
+
+ WRITE(numout,*) 'rst_read : modify restart to fit new geometry'
+ WRITE(numout,*) '~~~~~~~~'
+
+ CALL iom_get( numror, jpdom_autoglo, 'tmask' , ztmask_b ) ! need to extrapolate T/S
+ CALL iom_get( numror, jpdom_autoglo, 'umask' , zumask_b ) ! need to correct barotropic velocity
+ CALL iom_get( numror, jpdom_autoglo, 'vmask' , zvmask_b ) ! need to correct barotropic velocity
+ CALL iom_get( numror, jpdom_autoglo, 'smask' , zsmask_b ) ! need to correct barotropic velocity
+ CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', ze3t_b(:,:,:) ) ! need to compute temperature correction
+ CALL iom_get( numror, jpdom_autoglo, 'fse3u_n', ze3u_b(:,:,:) ) ! need to compute volume correction ????
+ CALL iom_get( numror, jpdom_autoglo, 'fse3v_n', ze3v_b(:,:,:) ) ! need to compute volume correction ????
+ CALL iom_get( numror, jpdom_autoglo, 'fsdepw_n', zdepw_b(:,:,:) ) ! need to compute volume correction ????
+ !! ! Extrapolation/interpolation of modify cell and new cells ... (maybe do it later after domvvl)
+ CALL rst_iscpl_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b, htsc_iscpl, hdiv_iscpl )
+ IF( nmsh /= 0 .AND. ln_iscpl ) CALL dom_wri ! Create a domain file
+
+ cfile='correction'
+ cfile = TRIM( cfile )
+ CALL iom_open ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )
+ CALL iom_rstput( 0, 0, inum0, 'vol_cor', hdiv_iscpl(:,:,:) )
+ CALL iom_rstput( 0, 0, inum0, 'tem_cor', htsc_iscpl(:,:,:,jp_tem) )
+ CALL iom_rstput( 0, 0, inum0, 'sal_cor', htsc_iscpl(:,:,:,jp_sal) )
+ CALL iom_rstput( 0, 0, inum0, 'e3t_diff', (e3t_0(:,:,:)*tmask(:,:,:) - ze3t_b(:,:,:)*ztmask_b(:,:,:))*tmask(:,:,:) )
+ CALL iom_close ( inum0 )
+
+ CALL wrk_dealloc(jpi,jpj,jpk, ztmask_b,zumask_b,zvmask_b ) ! no need of this variable afterwards (only used for interpolation)
+ CALL wrk_dealloc(jpi,jpj,jpk, ze3t_b ,ze3u_b ,ze3v_b ) ! no need of this variable afterwards (only used for interpolation)
+ CALL wrk_dealloc(jpi,jpj,jpk, zdepw_b )
+ CALL wrk_dealloc(jpi,jpj, zsmask_b )
+ END IF
+ IF( neuler == 0 ) THEN ! Euler restart (neuler=0)
+ tsb (:,:,:,:) = tsn (:,:,:,:) ! all before fields set to now values
+ ub (:,:,:) = un (:,:,:)
+ vb (:,:,:) = vn (:,:,:)
+ rotb (:,:,:) = rotn (:,:,:)
+ hdivb(:,:,:) = hdivn(:,:,:)
+ sshb (:,:) = sshn (:,:)
+
+ IF( lk_vvl ) THEN
+ DO jk = 1, jpk
+ fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
+ END DO
+ ENDIF
+
+ IF( lk_lim3 .AND. .NOT. lk_vvl ) THEN
+ DO jk = 1, jpk
+ fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
+ END DO
+ ENDIF
+
+ ENDIF
+ !
+ IF( lk_lim3 ) THEN
+ CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev' , fraqsr_1lev )
+ ENDIF
+
+ END SUBROUTINE rst_iscpl
+
+ SUBROUTINE rst_iscpl_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b, pts_flx, pvol_flx)
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE rst_iscpl ***
+ !!
+ !! ** Purpose : compute new tn, sn, un, vn and sshn in case of evolving geometry of ice shelves
+ !! compute 2d fields of heat, salt and volume correction
+ !!
+ !! ** Method : tn, sn : extrapolation from neigbourg cells
+ !! un, vn : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
+ !!----------------------------------------------------------------------
+ !!######## ~ cp/paste pyton script for this routine
+ REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b !! mask before
+ REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: pe3t_b , pe3u_b , pe3v_b !! scale factor before
+ REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: pdepw_b
+ REAL(wp), DIMENSION(:,: ), INTENT(in ) :: psmask_b !! mask before
+ REAL(wp), DIMENSION(:,:,:,:), INTENT(out) :: pts_flx !! corrective flux to have tracer conservation
+ REAL(wp), DIMENSION(:,:,: ), INTENT(out) :: pvol_flx !! corrective flux to have volume conservation
+ !!
+ INTEGER :: ji, jj, ii, ij, jk, iz !! loop index
+ INTEGER :: iip1, iim1, ijp1, ijm1, jkp1, jkm1
+ INTEGER :: ierror, inum ! temporary integer
+ INTEGER :: ios ! Local integer output status for namelist read
+ !!
+ REAL(wp):: summsk, r1_tiscpl, zsum, zsum1, zarea, zsumn, zsumb
+ REAL(wp):: ziip1_ratio, ziim1_ratio, zijp1_ratio, zijm1_ratio
+ REAL(wp):: zdz, zdzm1, zdzp1
+ !!
+ REAL(wp), DIMENSION(:,: ), POINTER :: zdmask , zdsmask, zvcorr, zucorr, zde3t
+ REAL(wp), DIMENSION(:,: ), POINTER :: zbub , zbvb , zbun , zbvn
+ REAL(wp), DIMENSION(:,: ), POINTER :: zssh0 , zssh1, hu1, hv1
+ REAL(wp), DIMENSION(:,: ), POINTER :: zsmask0, zsmask1
+ REAL(wp), DIMENSION(:,:,: ), POINTER :: ztmask0, ztmask1, ztrp
+ REAL(wp), DIMENSION(:,:,: ), POINTER :: zwmaskn, zwmaskb, ztmp3d
+ REAL(wp), DIMENSION(:,:,:,:), POINTER :: zts0
+ !
+ REAL(wp), DIMENSION(: ), ALLOCATABLE :: zlon, zlat
+ REAL(wp), DIMENSION(: ), ALLOCATABLE :: zcorr_vol, zcorr_tem, zcorr_sal
+ INTEGER , DIMENSION(: ), ALLOCATABLE :: ixpts, iypts, izpts, vnpts
+ INTEGER :: jpts, npts
+ !
+ NAMELIST/namsbc_iscpl/rn_fiscpl,ln_hfb
+ !!----------------------------------------------------------------------
+ ! ! ============
+ ! ! Namelist
+ ! ! ============
+ !
+ rn_fiscpl = 2480
+ ln_hfb = .FALSE.
+ REWIND( numnam_ref ) ! Namelist namsbc_iscpl in reference namelist : Ice sheet coupling
+ READ ( numnam_ref, namsbc_iscpl, IOSTAT = ios, ERR = 901)
+901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_iscpl in reference namelist', lwp )
+
+ REWIND( numnam_cfg ) ! Namelist namsbc_iscpl in configuration namelist : Ice Sheet coupling
+ READ ( numnam_cfg, namsbc_iscpl, IOSTAT = ios, ERR = 902 )
+902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_iscpl in configuration namelist', lwp )
+ IF(lwm) WRITE ( numond, namsbc_iscpl )
+ !
+ PRINT *, 'rn_iscpl .. =',rn_fiscpl, ln_hfb
+
+ CALL wrk_alloc(jpi,jpj,jpk,2, zts0 )
+ CALL wrk_alloc(jpi,jpj,jpk, ztmask0, ztmask1 , ztrp, ztmp3d )
+ CALL wrk_alloc(jpi,jpj,jpk, zwmaskn, zwmaskb )
+ CALL wrk_alloc(jpi,jpj, zsmask0, zsmask1 )
+ CALL wrk_alloc(jpi,jpj, zdmask , zdsmask, zvcorr, zucorr, zde3t)
+ CALL wrk_alloc(jpi,jpj, zbub , zbvb , zbun , zbvn )
+ CALL wrk_alloc(jpi,jpj, zssh0 , zssh1, hu1, hv1 )
+! mask value to be sure
+ tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:)
+ tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:)
+!=============================================================================
+
+ ! compute wmask
+ zwmaskn(:,:,1) = tmask (:,:,1)
+ zwmaskb(:,:,1) = ptmask_b(:,:,1)
+ DO jk = 2,jpk
+ zwmaskn(:,:,jk) = tmask (:,:,jk) * tmask (:,:,jk-1)
+ zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
+ END DO
+
+
+!=============================================================================
+
+ ! compute ssh
+ zsmask0(:,:) = psmask_b(:,:)
+ zsmask1(:,:) = psmask_b(:,:)
+
+ ! compute new ssh if we open a full water column (mean of the closest neigbourgs)
+ sshb (:,:)=sshn(:,:)
+ zssh0(:,:)=sshn(:,:)
+ DO iz = 1,10 ! need to be tuned (configuration dependent)
+ zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
+ DO ii = 2,jpi-1
+ DO ij = 2,jpj-1
+ iip1=ii+1; iim1=ii-1;
+ ijp1=ij+1; ijm1=ij-1;
+ summsk=(zsmask0(iip1,ij)+zsmask0(iim1,ij)+zsmask0(ii,ijp1)+zsmask0(ii,ijm1))
+ IF (zdsmask(ii,ij)==1 .AND. summsk .NE. 0) THEN
+ PRINT *, 'add ssh to 1 cell',ii,ij,narea
+ sshn(ii,ij)=(zssh0(iip1,ij)*zsmask0(iip1,ij) &
+ & +zssh0(iim1,ij)*zsmask0(iim1,ij) &
+ & +zssh0(ii,ijp1)*zsmask0(ii,ijp1) &
+ & +zssh0(ii,ijm1)*zsmask0(ii,ijm1))/summsk
+ zsmask1(ii,ij)=1
+ END IF
+ END DO
+ END DO
+ CALL lbc_lnk(sshn,'T',1.)
+ CALL lbc_lnk(zsmask1,'T',1.)
+ zssh0 = sshn
+ zsmask0 = zsmask1
+ END DO
+ sshn(:,:) = sshn(:,:) * ssmask(:,:)
+
+!=============================================================================
+! WARNING we cannot used glob_sum for before time variable (because glob_sum use tmask_i). Need to mask the halo and glob_sum_full
+ ztmp3d(:,:,:) = 0.0
+ ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
+ & - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
+ zsum = glob_sum_full(ztmp3d)
+ IF (lwp) PRINT *, 'total volume correction 00 = ',zsum
+ IF ( lk_vvl ) THEN
+! compute fse3t_n
+ DO jk = 1,jpk
+ fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1 + sshn(:,:) / ( ht_0(:,:) + 1. - ssmask(:,:) ) * tmask(:,:,jk) )
+ END DO
+ ztmp3d(:,:,:) = 0.0
+ ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
+ & - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
+ zsum = glob_sum_full(ztmp3d)
+ IF (lwp) PRINT *, 'total volume correction 01 = ',zsum
+
+! compute fse3u/v ... (call interpolation vvl)
+ ! Reconstruction of all vertical scale factors at now and before time steps
+ ! =============================================================================
+ ! Horizontal scale factor interpolations
+ ! --------------------------------------
+ CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
+ CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
+ CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
+ ! Vertical scale factor interpolations
+ ! ------------------------------------
+ CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' )
+ CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
+ CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
+ ! t- and w- points depth
+ ! ----------------------
+ fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
+ fsdepw_n(:,:,1) = 0.0_wp
+ fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
+ DO jj = 1,jpj
+ DO ji = 1,jpi
+ DO jk = 2,mikt(ji,jj)-1
+ fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
+ fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
+ fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
+ END DO
+ IF (mikt(ji,jj) .GT. 1) THEN
+ jk = mikt(ji,jj)
+ fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk)
+ fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
+ fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk ) - sshn (ji,jj)
+ END IF
+ DO jk = mikt(ji,jj)+1, jpk
+ fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)
+ fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
+ fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk ) - sshn (ji,jj)
+ END DO
+ END DO
+ END DO
+
+ hu(:,:) = 0._wp ! Ocean depth at U-points
+ hv(:,:) = 0._wp ! Ocean depth at V-points
+ ht(:,:) = 0._wp ! Ocean depth at T-points
+ DO jk = 1, jpkm1
+ hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
+ hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
+ ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+ END DO
+ ! ! Inverse of the local depth
+ hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
+ hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
+
+ END IF
+ ztmp3d(:,:,:) = 0.0
+ ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
+ & - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
+ zsum = glob_sum_full(ztmp3d)
+ IF (lwp) PRINT *, 'total volume correction 02 = ',zsum
+
+!=============================================================================
+
+! compute velocity
+! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
+ ub(:,:,:)=un(:,:,:)
+ vb(:,:,:)=vn(:,:,:)
+ DO jk = 1,jpk
+ DO ii = 1,jpi
+ DO ij = 1,jpj
+ un(ii,ij,jk) = ub(ii,ij,jk)*pe3u_b(ii,ij,jk)*pumask_b(ii,ij,jk)/fse3u_n(ii,ij,jk)*umask(ii,ij,jk);
+ vn(ii,ij,jk) = vb(ii,ij,jk)*pe3v_b(ii,ij,jk)*pvmask_b(ii,ij,jk)/fse3v_n(ii,ij,jk)*vmask(ii,ij,jk);
+ !IF ( umask(ii,ij,1) == 1 .AND. pumask(ii,ij,1) == 0 ) un(ii,ij,jk) = ub(ii,ij,jk)
+ !IF ( vmask(ii,ij,1) == 1 .AND. pvmask(ii,ij,1) == 0 ) vn(ii,ij,jk) = vb(ii,ij,jk)
+ END DO
+ END DO
+ END DO
+! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
+ ! compute barotropic velocity now and after
+ ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:);
+ zbub(:,:) = SUM(ztrp,DIM=3)
+ ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:);
+ zbvb(:,:) = SUM(ztrp,DIM=3)
+ ztrp(:,:,:) = un(:,:,:)*fse3u_n(:,:,:);
+ zbun(:,:) = SUM(ztrp,DIM=3)
+ ztrp(:,:,:) = vn(:,:,:)*fse3v_n(:,:,:);
+ zbvn(:,:) = SUM(ztrp,DIM=3)
+! Already know ????????
+ hu1=0.0_wp ;
+ hv1=0.0_wp ;
+ DO jk = 1,jpk
+ hu1(:,:) = hu1(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
+ hv1(:,:) = hv1(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
+ END DO
+ zucorr = 0._wp
+ zvcorr = 0._wp
+ DO ii = 1,jpi
+ DO ij = 1,jpj
+ IF (zbun(ii,ij) .NE. zbub(ii,ij) .AND. hu1(ii,ij) .NE. 0._wp ) THEN
+ zucorr(ii,ij) = (zbun(ii,ij) - zbub(ii,ij))/hu1(ii,ij)
+ END IF
+ IF (zbvn(ii,ij) .NE. zbvb(ii,ij) .AND. hv1(ii,ij) .NE. 0._wp ) THEN
+ zvcorr(ii,ij) = (zbvn(ii,ij) - zbvb(ii,ij))/hv1(ii,ij)
+ END IF
+ END DO
+ END DO
+ DO jk = 1,jpk
+ un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
+ vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)
+ END DO
+
+!=============================================================================
+ ! compute temp and salt
+ ! compute new tn and sn if we open a new cell
+ tsb (:,:,:,:) = tsn(:,:,:,:)
+ zts0(:,:,:,:) = tsn(:,:,:,:)
+ ztmask1(:,:,:) = ptmask_b(:,:,:)
+ ztmask0(:,:,:) = ptmask_b(:,:,:)
+ DO iz = 1,10
+ DO jk = 1,jpk-1
+ zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
+ DO ii = 2,jpi-1
+ DO ij = 2,jpj-1
+ iip1=ii+1; iim1=ii-1;
+ ijp1=ij+1; ijm1=ij-1;
+ summsk=(ztmask0(iip1,ij,jk)+ztmask0(iim1,ij,jk)+ztmask0(ii,ijp1,jk)+ztmask0(ii,ijm1,jk))
+ !! horizontal extrapolation
+ IF (zdmask(ii,ij)==1 .AND. summsk .NE. 0) THEN
+ tsn(ii,ij,jk,1)=( zts0(iip1,ij,jk,1)*ztmask0(iip1,ij,jk) &
+ & +zts0(iim1,ij,jk,1)*ztmask0(iim1,ij,jk) &
+ & +zts0(ii,ijp1,jk,1)*ztmask0(ii,ijp1,jk) &
+ & +zts0(ii,ijm1,jk,1)*ztmask0(ii,ijm1,jk))/summsk
+ tsn(ii,ij,jk,2)=( zts0(iip1,ij,jk,2)*ztmask0(iip1,ij,jk) &
+ & +zts0(iim1,ij,jk,2)*ztmask0(iim1,ij,jk) &
+ & +zts0(ii,ijp1,jk,2)*ztmask0(ii,ijp1,jk) &
+ & +zts0(ii,ijm1,jk,2)*ztmask0(ii,ijm1,jk))/summsk
+ ztmask1(ii,ij,jk)=1
+ END IF
+ !! vertical extrapolation if horizontal extra failed
+ IF (zdmask(ii,ij)==1 .AND. summsk==0) THEN
+ jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
+ summsk=(ztmask0(ii,ij,jkm1)+ztmask0(ii,ij,jkp1))
+ IF (zdmask(ii,ij)==1 .AND. summsk .NE. 0 ) THEN
+ tsn(ii,ij,jk,1)=( zts0(ii,ij,jkp1,1)*ztmask0(ii,ij,jkp1) &
+ & +zts0(ii,ij,jkm1,1)*ztmask0(ii,ij,jkm1))/summsk
+ tsn(ii,ij,jk,2)=( zts0(ii,ij,jkp1,2)*ztmask0(ii,ij,jkp1) &
+ & +zts0(ii,ij,jkm1,2)*ztmask0(ii,ij,jkm1))/summsk
+ ztmask1(ii,ij,jk)=1
+ END IF
+ END IF
+ !! horizontal corner extrapolation if horizontal and vertical failed
+ IF (zdmask(ii,ij)==1 .AND. summsk==0) THEN
+ iip1=ii+1; iim1=ii-1;
+ ijp1=ij+1; ijm1=ij-1;
+ summsk=(ztmask0(iip1,ijp1,jk)+ztmask0(iim1,ijm1,jk)+ztmask0(iim1,ijp1,jk)+ztmask0(iip1,ijm1,jk))
+ IF (zdmask(ii,ij)==1 .AND. summsk .NE. 0) THEN
+ tsn(ii,ij,jk,1)=( zts0(iip1,ijp1,jk,1)*ztmask0(iip1,ijp1,jk) &
+ & +zts0(iim1,ijm1,jk,1)*ztmask0(iim1,ijm1,jk) &
+ & +zts0(iim1,ijp1,jk,1)*ztmask0(iim1,ijp1,jk) &
+ & +zts0(iip1,ijm1,jk,1)*ztmask0(iip1,ijm1,jk))/summsk
+ tsn(ii,ij,jk,2)=( zts0(iip1,ijp1,jk,2)*ztmask0(iip1,ijp1,jk) &
+ & +zts0(iim1,ijm1,jk,2)*ztmask0(iim1,ijm1,jk) &
+ & +zts0(iim1,ijp1,jk,2)*ztmask0(iim1,ijp1,jk) &
+ & +zts0(iip1,ijm1,jk,2)*ztmask0(iip1,ijm1,jk))/summsk
+ ztmask1(ii,ij,jk)=1
+ END IF
+ END IF
+ END DO
+ END DO
+ END DO
+ CALL lbc_lnk(tsn(:,:,:,1),'T',1.)
+ CALL lbc_lnk(tsn(:,:,:,2),'T',1.)
+ CALL lbc_lnk(ztmask1,'T',1.)
+ zts0(:,:,:,:) = tsn(:,:,:,:)
+ ztmask0 = ztmask1
+ END DO
+ tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
+ tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
+
+ ! CHECK heat
+ zsumn = glob_sum( tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) * zwmaskb(:,:,:) * zwmaskn(:,:,:))
+ zsumb = glob_sum( tsb(:,:,:,jp_tem) * pe3t_b(:,:,:) * zwmaskb(:,:,:) * zwmaskn(:,:,:))
+ IF (lwp) PRINT *, 'CHECK tsn = ',zsumn, zsumb
+
+ ! compute new T/S (interpolation) because of vvl
+ IF ( lk_vvl .AND. .FALSE. ) THEN
+ !IF ( lk_vvl ) THEN
+ DO jk = 2,jpk-1
+ DO jj = 1,jpj
+ DO ji = 1,jpi
+ IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1.0_wp ) THEN
+ !compute weight
+ zdzp1 = MAX(0._wp,fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
+ zdz = fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk )
+ zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk ) - fsdepw_n(ji,jj,jk ))
+ IF (zdz .LT. 0.0_wp) THEN
+ PRINT *, 'ERROR dz n ', ji,jj,jk,zdz,fsdepw_n(ji,jj,jk+1),fsdepw_n(ji,jj,jk),fsdepw_n(ji,jj,jk-1)
+ PRINT *, 'ERROR dz n = ',fse3t_n (ji,jj,jk+1),fse3t_n (ji,jj,jk),fse3t_n (ji,jj,jk-1), sshn(ji,jj)
+ PRINT *, 'ERROR dz b = ',pdepw_b(ji,jj,jk+1),pdepw_b(ji,jj,jk),pdepw_b(ji,jj,jk-1)
+ PRINT *, 'ERROR dz b = ',fse3t_b (ji,jj,jk+1),fse3t_b (ji,jj,jk),fse3t_b (ji,jj,jk-1), sshb(ji,jj)
+ PRINT *, 'ERROR dz 0 = ', e3t_0 (ji,jj,jk+1), e3t_0 (ji,jj,jk), e3t_0 (ji,jj,jk-1)
+ PRINT *, 'ERROR dz n = ', tmask (ji,jj,jk+1), tmask (ji,jj,jk), tmask (ji,jj,jk-1)
+ PRINT *, 'ERROR dz n = ', zwmaskn(ji,jj,jk+1), zwmaskn(ji,jj,jk), zwmaskn(ji,jj,jk-1)
+ PRINT *, 'ERROR dz b = ', ptmask_b(ji,jj,jk+1), ptmask_b(ji,jj,jk), ptmask_b(ji,jj,jk-1)
+ PRINT *, 'ERROR dz b = ', zwmaskb(ji,jj,jk+1), zwmaskb(ji,jj,jk), zwmaskb(ji,jj,jk-1)
+ PRINT *, 'ERROR dz b = ', gdepw_0(ji,jj,jk+1), gdepw_0(ji,jj,jk), gdepw_0(ji,jj,jk-1)
+ STOP
+ END IF
+ tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem) &
+ & + zdz *tsb(ji,jj,jk ,jp_tem) &
+ & + zdzm1*tsb(ji,jj,jk-1,jp_tem) )/fse3t_n(ji,jj,jk)
+ tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal) &
+ & + zdz *tsb(ji,jj,jk ,jp_sal) &
+ & + zdzm1*tsb(ji,jj,jk-1,jp_sal) )/fse3t_n(ji,jj,jk)
+ END IF
+ END DO
+ END DO
+ END DO
+ END IF
+
+ ! CHECK heat
+ ztmp3d(:,:,:) = 0.0
+ ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
+ & - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
+ zsum = glob_sum_full(ztmp3d)
+ IF (lwp) PRINT *, 'total volume correction 03 = ',zsum
+
+ WHERE (tmask(:,:,:) == 1.0 .AND. tsn(:,:,:,2) == 0._wp)
+ tsn(:,:,:,2)= -99._wp
+ tmask(:,:,:) = 0.0
+ umask(:,:,:) = 0.0
+ vmask(:,:,:) = 0.0
+ END WHERE
+
+ WHERE (SUM(tmask,dim=3) == 0)
+ mbkt(:,:)=1 ; mbku(:,:)=1 ; mbkv(:,:)=1
+ mikt(:,:)=1 ; miku(:,:)=1 ; mikv(:,:)=1
+ END WHERE
+ ! compute new tn and sn if we close cell
+! nothing to do
+
+!=============================================================================
+! put this stuff in a subroutine ????
+ !IF ( ln_hfb ) CALL scb_iscpl_cons(pvol_flx, pts_flx, pe3t_b, zssh0)
+ ! initialisation
+ zde3t (:,:) = 0.0_wp
+ pvol_flx(:,:,: ) = 0.0_wp
+ pts_flx (:,:,:,:) = 0.0_wp
+ IF ( ln_hfb ) THEN
+ zsum = glob_sum_full(pvol_flx(:,:,:) ) * rn_fiscpl * rn_rdt
+ IF (lwp) PRINT *, 'total volume correction 0 = ',zsum
+ zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
+ IF (lwp) PRINT *, 'total heat correction 0 = ',zsum
+ zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
+ IF (lwp) PRINT *, 'total salt correction 0 = ',zsum
+
+ ! mask tsn and tsb (should be useless)
+ tsb(:,:,:,jp_tem)=tsb(:,:,:,jp_tem)*ptmask_b(:,:,:); tsn(:,:,:,jp_tem)=tsn(:,:,:,jp_tem)*tmask(:,:,:);
+ tsb(:,:,:,jp_sal)=tsb(:,:,:,jp_sal)*ptmask_b(:,:,:); tsn(:,:,:,jp_sal)=tsn(:,:,:,jp_sal)*tmask(:,:,:);
+ ! why tsn = tsb ????? (bug ???????)
+ ! diagnose non conservation of heat, salt and volume
+ r1_tiscpl = 1._wp / (rn_fiscpl * rn_rdt)
+ zssh0(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
+ IF ( lk_vvl ) zssh0 = 0.0_wp
+ DO jk = 1,jpk-1
+ DO ii = 2,jpi-1
+ DO ij = 2,jpj-1
+ ! volume differences
+ zde3t(ii,ij) = fse3t_n(ii,ij,jk) * tmask(ii,ij,jk) - pe3t_b(ii,ij,jk) * ptmask_b(ii,ij,jk);
+
+ ! shh changes
+ IF ( ptmask_b(ii,ij,jk) == 1 .OR. tmask(ii,ij,jk) == 1 ) THEN
+ zde3t(ii,ij) = zde3t(ii,ij) + zssh0(ii,ij)
+ zssh0(ii,ij) = 0._wp
+ END IF
+
+ ! ocean cell now
+! IF (tmask(ii,ij,jk) == 1._wp) THEN
+ ! case where we open, enlarge or thin a cell :
+ pvol_flx(ii,ij,jk) = zde3t(ii,ij) * r1_tiscpl
+ pts_flx (ii,ij,jk,jp_sal)= tsn(ii,ij,jk,jp_sal) * zde3t(ii,ij) * r1_tiscpl
+ pts_flx (ii,ij,jk,jp_tem)= tsn(ii,ij,jk,jp_tem) * zde3t(ii,ij) * r1_tiscpl
+! END IF
+ END DO
+ END DO
+ END DO
+ ! glob_sum_full because with glob summ some data can be masked. WARNING the halo have to be set at 0
+ PRINT *, 'test ', narea, SUM(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt, SUM(pvol_flx(2:jpi-1,2:jpj-1,:)) * rn_fiscpl * rn_rdt
+ zsum = glob_sum_full(pvol_flx(:,:,:) ) * rn_fiscpl * rn_rdt
+ IF (lwp) PRINT *, 'total volume correction 1 = ',zsum
+ zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
+ IF (lwp) PRINT *, 'total heat correction 1 = ',zsum
+ zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
+ IF (lwp) PRINT *, 'total salt correction 1 = ',zsum
+
+ zssh0(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
+ IF ( lk_vvl ) zssh0 = 0.0_wp
+ DO jk = 1,jpk-1
+ DO ii = 2,jpi-1
+ DO ij = 2,jpj-1
+ ! volume differences
+ zde3t(ii,ij) = fse3t_n(ii,ij,jk) * tmask(ii,ij,jk) - pe3t_b(ii,ij,jk) * ptmask_b(ii,ij,jk);
+
+ ! shh changes
+ !IF ( ( ptmask_b(ii,ij,jk) == 1 .OR. tmask(ii,ij,jk) == 1 ) .AND. zssh0(ii,ij) .NE. 0._wp) THEN
+ IF ( ptmask_b(ii,ij,jk) == 1 .OR. tmask(ii,ij,jk) == 1 ) THEN
+ zde3t(ii,ij) = zde3t(ii,ij) + zssh0(ii,ij)
+ zssh0(ii,ij) = 0._wp
+ END IF
+
+ ! ocean cell before and mask cell now
+ IF ( tmask(ii,ij,jk) == 0._wp .AND. ptmask_b(ii,ij,jk) == 1._wp ) THEN
+ ! case where we close a cell and adjacent cell open
+ pvol_flx(ii,ij,jk) = zde3t(ii,ij) * r1_tiscpl
+ pts_flx (ii,ij,jk,jp_sal)= tsb(ii,ij,jk,jp_sal) * zde3t(ii,ij) * r1_tiscpl
+ pts_flx (ii,ij,jk,jp_tem)= tsb(ii,ij,jk,jp_tem) * zde3t(ii,ij) * r1_tiscpl
+
+ iip1=ii+1 ; iim1=ii-1 ; ijp1=ij+1 ; ijm1=ij-1 ;
+
+ zsum = e12t(ii ,ijp1) * tmask(ii ,ijp1,jk) + e12t(ii ,ijm1) * tmask(ii ,ijm1,jk) &
+ & + e12t(iim1,ij ) * tmask(iim1,ij ,jk) + e12t(iip1,ij ) * tmask(iip1,ij ,jk)
+
+ IF ( zsum .NE. 0._wp ) THEN
+ ziip1_ratio = e12t(iip1,ij ) * tmask(iip1,ij ,jk) / zsum
+ ziim1_ratio = e12t(iim1,ij ) * tmask(iim1,ij ,jk) / zsum
+ zijp1_ratio = e12t(ii ,ijp1) * tmask(ii ,ijp1,jk) / zsum
+ zijm1_ratio = e12t(ii ,ijm1) * tmask(ii ,ijm1,jk) / zsum
+
+ pvol_flx(ii ,ijp1,jk ) = pvol_flx(ii ,ijp1,jk ) + pvol_flx(ii,ij,jk ) * zijp1_ratio
+ pvol_flx(ii ,ijm1,jk ) = pvol_flx(ii ,ijm1,jk ) + pvol_flx(ii,ij,jk ) * zijm1_ratio
+ pvol_flx(iip1,ij ,jk ) = pvol_flx(iip1,ij ,jk ) + pvol_flx(ii,ij,jk ) * ziip1_ratio
+ pvol_flx(iim1,ij ,jk ) = pvol_flx(iim1,ij ,jk ) + pvol_flx(ii,ij,jk ) * ziim1_ratio
+ pts_flx (ii ,ijp1,jk,jp_sal) = pts_flx (ii ,ijp1,jk,jp_sal) + pts_flx (ii,ij,jk,jp_sal) * zijp1_ratio
+ pts_flx (ii ,ijm1,jk,jp_sal) = pts_flx (ii ,ijm1,jk,jp_sal) + pts_flx (ii,ij,jk,jp_sal) * zijm1_ratio
+ pts_flx (iip1,ij ,jk,jp_sal) = pts_flx (iip1,ij ,jk,jp_sal) + pts_flx (ii,ij,jk,jp_sal) * ziip1_ratio
+ pts_flx (iim1,ij ,jk,jp_sal) = pts_flx (iim1,ij ,jk,jp_sal) + pts_flx (ii,ij,jk,jp_sal) * ziim1_ratio
+ pts_flx (ii ,ijp1,jk,jp_tem) = pts_flx (ii ,ijp1,jk,jp_tem) + pts_flx (ii,ij,jk,jp_tem) * zijp1_ratio
+ pts_flx (ii ,ijm1,jk,jp_tem) = pts_flx (ii ,ijm1,jk,jp_tem) + pts_flx (ii,ij,jk,jp_tem) * zijm1_ratio
+ pts_flx (iip1,ij ,jk,jp_tem) = pts_flx (iip1,ij ,jk,jp_tem) + pts_flx (ii,ij,jk,jp_tem) * ziip1_ratio
+ pts_flx (iim1,ij ,jk,jp_tem) = pts_flx (iim1,ij ,jk,jp_tem) + pts_flx (ii,ij,jk,jp_tem) * ziim1_ratio
+
+ ! set to 0 the cell we distributed over neigbourg cells
+ pvol_flx(ii,ij,jk ) = 0._wp
+ pts_flx (ii,ij,jk,jp_sal) = 0._wp
+ pts_flx (ii,ij,jk,jp_tem) = 0._wp
+ ELSE IF (zsum .EQ. 0._wp ) THEN
+ ! case where we close a cell and no adjacent cell open
+ IF ( tmask(ii,ij,jk+1) .EQ. 1._wp ) THEN
+ pvol_flx(ii,ij,jk+1) = pvol_flx(ii,ij,jk+1) + pvol_flx(ii,ij,jk)
+ pts_flx (ii,ij,jk+1,jp_sal)= pts_flx (ii,ij,jk+1,jp_sal) + pts_flx (ii,ij,jk,jp_sal)
+ pts_flx (ii,ij,jk+1,jp_tem)= pts_flx (ii,ij,jk+1,jp_tem) + pts_flx (ii,ij,jk,jp_tem)
+
+ ! set to 0 the cell we distributed over neigbourg cells
+ pvol_flx(ii,ij,jk ) = 0._wp
+ pts_flx (ii,ij,jk,jp_sal) = 0._wp
+ pts_flx (ii,ij,jk,jp_tem) = 0._wp
+ ELSE
+ ! case no adjacent cell on the horizontal and on the vertical
+ PRINT *, 'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal'
+ PRINT *, ' ',mig(ii),' ',mjg(ij),' ',jk
+ PRINT *, ' ',ii,' ',ij,' ',jk,' ',narea
+ PRINT *, ' conservation broken '
+! We deal with this point later.
+ END IF
+ END IF
+ END IF
+ END DO
+ END DO
+ END DO
+
+ zsum = glob_sum_full(pvol_flx(:,:,:) ) * rn_fiscpl * rn_rdt
+ IF (lwp) PRINT *, 'total volume correction 2 = ',zsum
+ zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
+ IF (lwp) PRINT *, 'total heat correction 2 = ',zsum
+ zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
+ IF (lwp) PRINT *, 'total salt correction 2 = ',zsum
+
+ ! allocation and initialisation of the list of problematic point
+ ALLOCATE(vnpts(jpnij))
+ vnpts(:)=0
+
+ ! fill narea location with the number of problematic point
+ DO jk = 1,jpk-1
+ DO ii = 2,jpi-1
+ DO ij = 2,jpj-1
+ IF ( ptmask_b(ii,ij,jk) == 1 .AND. SUM(tmask(ii-1:ii+1,ij,jk)) == 0 &
+ & .AND. SUM(tmask(ii,ij-1:ij+1,jk)) == 0 .AND. tmask(ii,ij,jk+1) == 0 ) THEN
+ vnpts(narea) = vnpts(narea) + 1
+ END IF
+ END DO
+ END DO
+ END DO
+
+ ! build array of total problematic point on each cpu (share to each cpu)
+ CALL mpp_max(vnpts,jpnij)
+
+ ! size of the new variable
+ npts = SUM(vnpts)
+
+ ! allocation of the coordiantes, correction, index vector for the problematic points
+ ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts))
+ ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20 ; zlat(:) = -1.0e20
+ zcorr_vol(:) = 0.0_wp
+ zcorr_sal(:) = 0.0_wp
+ zcorr_tem(:) = 0.0_wp
+
+ ! fill new variable
+ jpts = SUM(vnpts(1:narea-1))
+ DO jk = 1,jpk-1
+ DO ii = 2,jpi-1
+ DO ij = 2,jpj-1
+ IF ( ptmask_b(ii,ij,jk) == 1 .AND. SUM(tmask(ii-1:ii+1,ij,jk)) == 0 &
+ & .AND. SUM(tmask(ii,ij-1:ij+1,jk)) == 0 .AND. tmask(ii,ij,jk+1) == 0 ) THEN
+ jpts = jpts + 1
+ PRINT *, 'corrected point ', narea, ii, ij, jk, jpts
+ ixpts(jpts) = ii ; iypts(jpts) = ij ; izpts(jpts) = jk
+ zlon (jpts) = glamt(ii,ij) ; zlat (jpts) = gphit(ii,ij)
+ zcorr_vol(jpts) = pvol_flx(ii,ij,jk)
+ zcorr_sal(jpts) = pts_flx (ii,ij,jk,jp_sal)
+ zcorr_tem(jpts) = pts_flx (ii,ij,jk,jp_tem)
+ ! set flx to 0 (safer)
+ pvol_flx(ii,ij,jk ) = 0.0_wp
+ pts_flx (ii,ij,jk,jp_sal) = 0.0_wp
+ pts_flx (ii,ij,jk,jp_tem) = 0.0_wp
+ PRINT *, zcorr_vol(jpts)*rn_fiscpl*rn_rdt, zcorr_sal(jpts)*rn_fiscpl*rn_rdt, zcorr_tem(jpts)*rn_fiscpl*rn_rdt
+ END IF
+ END DO
+ END DO
+ END DO
+
+ ! build array of total problematic point on each cpu (share to each cpu)
+ CALL mpp_max(zlat ,npts)
+ CALL mpp_max(zlon ,npts)
+ CALL mpp_max(izpts,npts)
+
+ ! put correction term in the closest cell
+ PRINT *, 'corrected point1 ', narea, zlon, zlat, izpts
+ DO jpts = 1,npts
+ CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts))
+ PRINT *, 'corrected point2 ', narea, jpts, ixpts(jpts), iypts(jpts), izpts(jpts)
+ DO ij = mj0(iypts(jpts)),mj1(iypts(jpts))
+ DO ii = mi0(ixpts(jpts)),mi1(ixpts(jpts))
+ jk = izpts(jpts)
+ pvol_flx(ii,ij,jk) = pvol_flx(ii,ij,jk ) + zcorr_vol(jpts)
+ pts_flx (ii,ij,jk,jp_sal) = pts_flx (ii,ij,jk,jp_sal) + zcorr_sal(jpts)
+ pts_flx (ii,ij,jk,jp_tem) = pts_flx (ii,ij,jk,jp_tem) + zcorr_tem(jpts)
+ END DO
+ END DO
+ END DO
+ ! deallocate variable
+ DEALLOCATE(vnpts)
+ DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat)
+
+ ! add contribution store on the allo (lbclnk remove one of the contribution)
+ pvol_flx(:,:,: ) = pvol_flx(:,:,: ) * tmask(:,:,:)
+ pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:)
+ pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:)
+
+ CALL lbc_sum(pvol_flx(:,:,: ),'T',1.)
+ CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.)
+ CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.)
+
+ ! CHECK vol !!!!!!!!! warning tmask_i wrong if deals with before value, so glob_sum wrong for before value!!!!
+ zsumn = glob_sum ( fse3t_n(:,:,:) * tmask (:,:,:)) - glob_sum(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt
+ ztmp3d(:,:,:) = 0.0
+ ztmp3d(2:jpi-1,2:jpj-1,:) = pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
+ zsumb = glob_sum_full(ztmp3d)
+ zsum = glob_sum ( pvol_flx(:,:,:) * rn_fiscpl * rn_rdt)
+ IF (lwp) PRINT *, 'CHECK vol = ',zsumn, zsumb, zsumn - zsumb, zsum
+ ! CHECK salt
+ zsumn = glob_sum( tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) * tmask (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
+ zsumb = glob_sum( tsb(:,:,:,jp_sal) * pe3t_b(:,:,:) * ptmask_b(:,:,:))
+ zsum = glob_sum( pts_flx(:,:,:,jp_sal)*rn_fiscpl * rn_rdt)
+ IF (lwp) PRINT *, 'CHECK salt = ',zsumn, zsumb, zsumn - zsumb, zsum
+ ! CHECK heat
+ zsumn = glob_sum( tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) * tmask (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
+ zsumb = glob_sum( tsb(:,:,:,jp_tem) * pe3t_b(:,:,:) * ptmask_b(:,:,:))
+ zsum = glob_sum( pts_flx(:,:,:,jp_tem)*rn_fiscpl * rn_rdt)
+ IF (lwp) PRINT *, 'CHECK heat = ',zsumn, zsumb, zsumn - zsumb, zsum
+
+ END IF
+!=============================================================================
+ !! now T/S/U are compute in agreement with the new grid.
+ !! need to interpolate on vvl grid
+ !! next step is an euler time step
+ neuler = 0
+ !! set _b and _n variables equal
+ tsb (:,:,:,:)=tsn (:,:,:,:)
+ ub (:,:,: )=un (:,:,: )
+ vb (:,:,: )=vn (:,:,: )
+ sshb(:,: )=sshn(:,:)
+ !! set _b and _n vertical scale factor equal
+ fse3t_b (:,:,:)=fse3t_n (:,:,:)
+ fse3u_b (:,:,:)=fse3u_n (:,:,:)
+ fse3v_b (:,:,:)=fse3v_n (:,:,:)
+ IF ( lk_vvl ) THEN
+ fse3uw_b(:,:,:) = fse3uw_n(:,:,:)
+ fse3vw_b(:,:,:) = fse3vw_n(:,:,:)
+ !! set depth_b = depth_n
+ fsdept_b(:,:,:) = fsdept_n(:,:,:)
+ fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
+ !! set hu_b
+ hu_b(:,:) = hu(:,:)
+ hv_b(:,:) = hv(:,:)
+ hur_b(:,:) = hur(:,:)
+ hvr_b(:,:) = hvr(:,:)
+ END IF
+ !!
+ CALL wrk_dealloc(jpi,jpj,jpk,2, zts0 )
+ CALL wrk_dealloc(jpi,jpj,jpk, ztmask0, ztmask1 , ztrp )
+ CALL wrk_dealloc(jpi,jpj,jpk, zwmaskn, zwmaskb , ztmp3d )
+ CALL wrk_dealloc(jpi,jpj, zsmask0, zsmask1 )
+ CALL wrk_dealloc(jpi,jpj, zdmask , zdsmask, zvcorr, zucorr, zde3t)
+ CALL wrk_dealloc(jpi,jpj, zbub , zbvb , zbun , zbvn )
+ CALL wrk_dealloc(jpi,jpj, zssh0 , zssh1 , hu1 , hv1 )
+ END SUBROUTINE rst_iscpl_interpol
+
+ SUBROUTINE sbc_iscpl_div( phdivn )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE sbc_iscpl_div ***
+ !!
+ !! ** Purpose : update the horizontal divergence with the iscpl imbalance
+ !!
+ !! ** Method :
+ !! CAUTION : iscpl is positive (inflow) increase the
+ !! divergence and expressed in m/s
+ !!
+ !! ** Action : phdivn decreased by the runoff inflow
+ !!----------------------------------------------------------------------
+ REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: phdivn ! horizontal divergence
+ !!
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ !!----------------------------------------------------------------------
+ !
+ DO jk = 1, jpk
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / fse3t_n(ji,jj,jk)
+ END DO
+ END DO
+ END DO
+ !
+ END SUBROUTINE sbc_iscpl_div
+
+END MODULE sbc_iscpl
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 (revision 5619)
@@ -63,5 +63,5 @@
REAL(wp), PUBLIC, SAVE :: rcpi = 2000.0_wp ! phycst ?
- REAL(wp), PUBLIC, SAVE :: kappa = 1.54e-6_wp ! phycst ?
+ REAL(wp), PUBLIC, SAVE :: kappa = 0.0_wp ! phycst ?
REAL(wp), PUBLIC, SAVE :: rhoisf = 920.0_wp ! phycst ?
REAL(wp), PUBLIC, SAVE :: tsurf = -20.0_wp ! phycst ?
@@ -125,4 +125,6 @@
IF ( lwp ) WRITE(numout,*) ' ln_divisf = ', ln_divisf
IF ( lwp ) WRITE(numout,*) ' nn_gammablk = ', nn_gammablk
+ IF ( lwp ) WRITE(numout,*) ' rn_gammat0 = ', rn_gammat0
+ IF ( lwp ) WRITE(numout,*) ' rn_gammas0 = ', rn_gammas0
IF ( lwp ) WRITE(numout,*) ' rn_tfri2 = ', rn_tfri2
IF (ln_divisf) THEN ! keep it in the namelist ??? used true anyway as for runoff ? (PM)
@@ -150,7 +152,6 @@
!: read effective lenght (BG03)
IF (nn_isf == 2) THEN
- ! Read Data and save some integral values
+ cvarLeff = 'soLeff'
CALL iom_open( sn_Leff_isf%clname, inum )
- cvarLeff = 'soLeff' !: variable name for Efficient Length scale
CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
CALL iom_close(inum)
@@ -237,11 +238,11 @@
CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')
CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T')
- CALL sbc_isf_tbl(un(:,:,:),utbl(:,:),'U')
- CALL sbc_isf_tbl(vn(:,:,:),vtbl(:,:),'V')
+ CALL sbc_isf_tbl(un(:,:,:) ,utbl(:,:),'U')
+ CALL sbc_isf_tbl(vn(:,:,:) ,vtbl(:,:),'V')
! iom print
CALL iom_put('ttbl',ttbl(:,:))
CALL iom_put('stbl',stbl(:,:))
- CALL iom_put('utbl',utbl(:,:))
- CALL iom_put('vtbl',vtbl(:,:))
+ CALL iom_put('utbl',utbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:))
+ CALL iom_put('vtbl',vtbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:))
! compute fwf and heat flux
CALL sbc_isf_cav (kt)
@@ -296,6 +297,6 @@
!
! output
- CALL iom_put('qisf' , qisf)
- IF( iom_use('fwfisf') ) CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce )
+ IF( iom_use('qisf' ) ) CALL iom_put('qisf' , qisf)
+ IF( iom_use('fwfisf') ) CALL iom_put('fwfisf', fwfisf)
END IF
@@ -416,5 +417,5 @@
REAL(wp) :: zfwflx, zhtflx, zhtflx_b
REAL(wp) :: zgammat, zgammas
- REAL(wp) :: zeps = -1.e-20_wp !== Local constant initialization ==!
+ REAL(wp) :: zeps = 1.e-20_wp !== Local constant initialization ==!
INTEGER :: ji, jj ! dummy loop indices
INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
@@ -506,10 +507,10 @@
zeps1=rcp*rau0*zgammat
zeps2=lfusisf*rau0*zgammas
- zeps3=rhoisf*rcpi*kappa/risfdep(ji,jj)
+ zeps3=rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps)
zeps4=zlamb2+zlamb3*risfdep(ji,jj)
zeps6=zeps4-zti(ji,jj)
zeps7=zeps4-tsurf
zaqe=zlamb1 * (zeps1 + zeps3)
- zaqer=0.5/zaqe
+ zaqer=0.5/MIN(zaqe,-zeps)
zbqe=zeps1*zeps6+zeps3*zeps7-zeps2
zcqe=zeps2*stbl(ji,jj)
@@ -520,6 +521,6 @@
zfrz(ji,jj)=zeps4+zlamb1*zsfrz
-! zfwflx is upward water flux
- zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / zsfrz )
+! zfwflx is upward water flux (MAX usefull if kappa = 0
+ zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) )
! zhtflx is upward heat flux (out of ocean)
! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
@@ -527,5 +528,5 @@
! zwflx is upward water flux
! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
- zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj)
+!!!!!!!! zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj)
! test convergence and compute gammat
IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE.
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90 (revision 5619)
@@ -921,5 +921,5 @@
pn2(ji,jj,jk) = grav * ( zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) &
& - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) &
- & / fse3w(ji,jj,jk) * tmask(ji,jj,jk)
+ & / fse3w(ji,jj,jk) * wmask(ji,jj,jk)
END DO
END DO
@@ -1242,5 +1242,5 @@
!
rau0 = 1026._wp !: volumic mass of reference [kg/m3]
- rcp = 3991.86795711963_wp !: heat capacity [J/K]
+ rcp = 3974._wp !: heat capacity [J/K]
!
IF(lwp) THEN ! Control print
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 (revision 5619)
@@ -102,12 +102,13 @@
REAL(wp) :: zta, zsa ! local scalars
REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zts_dta
+ REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zts_dtadmp
!!----------------------------------------------------------------------
!
IF( nn_timing == 1 ) CALL timing_start( 'tra_dmp')
!
- CALL wrk_alloc( jpi, jpj, jpk, jpts, zts_dta )
- !
+ CALL wrk_alloc( jpi, jpj, jpk, jpts, zts_dta, zts_dtadmp )
! !== input T-S data at kt ==!
- CALL dta_tsd( kt, zts_dta ) ! read and interpolates T-S data at kt
+ CALL dta_tsd( kt, zts_dta, zts_dtadmp ) ! read and interpolates T-S data at kt
+ zts_dta=zts_dtadmp
!
SELECT CASE ( nn_zdmp ) !== type of damping ==!
@@ -175,5 +176,5 @@
& tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
!
- CALL wrk_dealloc( jpi, jpj, jpk, jpts, zts_dta )
+ CALL wrk_dealloc( jpi, jpj, jpk, jpts, zts_dta, zts_dtadmp )
!
IF( nn_timing == 1 ) CALL timing_stop( 'tra_dmp')
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90 (revision 5619)
@@ -22,4 +22,5 @@
USE sbcrnf ! River runoff
USE sbcisf ! Ice shelf
+ USE sbc_iscpl ! Ice sheet coupling
USE traqsr ! solar radiation penetration
USE trd_oce ! trends: ocean variables
@@ -117,8 +118,6 @@
INTEGER :: ji, jj, jk, jn ! dummy loop indices
INTEGER :: ikt, ikb
- INTEGER :: nk_isf
REAL(wp) :: zfact, z1_e3t, zdep
- REAL(wp) :: zalpha, zhk
- REAL(wp) :: zt_frz, zpress
+ REAL(wp) :: zt_frz, zpress
REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds
!!----------------------------------------------------------------------
@@ -220,5 +219,5 @@
!
IF( nn_isf > 0 ) THEN
- zfact = 0.5e0
+ zfact = 0.5_wp
DO jj = 2, jpj
DO ji = fs_2, fs_jpim1
@@ -233,5 +232,5 @@
! compute tfreez for the temperature correction (we add water at freezing temperature)
! zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04
- zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress )
+ zt_frz = -1.9_wp !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress )
! compute trend
tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) &
@@ -246,5 +245,5 @@
! compute tfreez for the temperature correction (we add water at freezing temperature)
! zpress = grav*rau0*fsdept(ji,jj,ikb)*1.e-04
- zt_frz = -1.9 !eos_fzp( tsn(ji,jj,ikb,jp_sal), zpress )
+ zt_frz = -1.9_wp !eos_fzp( tsn(ji,jj,ikb,jp_sal), zpress )
! compute trend
tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem) &
@@ -288,5 +287,23 @@
ENDIF
- IF( l_trdtra ) THEN ! send trends for further diagnostics
+ !----------------------------------------
+ ! Ice Sheet coupling imbalance correction to have conservation
+ !----------------------------------------
+ !
+ IF( ln_iscpl .AND. ln_hfb) THEN ! input of heat and salt due to river runoff
+ DO jk = 1,jpk
+ DO jj = 2, jpj
+ DO ji = fs_2, fs_jpim1
+ zdep = 1._wp / fse3t_n(ji,jj,jk)
+ tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) &
+ & * zdep
+ tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) &
+ & * zdep
+ END DO
+ END DO
+ END DO
+ ENDIF
+
+ IF( l_trdtra ) THEN ! save the horizontal diffusive trends for further diagnostics
ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
Index: /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90
===================================================================
--- /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90 (revision 5618)
+++ /branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90 (revision 5619)
@@ -25,4 +25,5 @@
PUBLIC glob_sum ! used in many places
+ PUBLIC glob_sum_full ! used in many places
PUBLIC DDPDD ! also used in closea module
PUBLIC glob_min, glob_max
@@ -34,4 +35,7 @@
MODULE PROCEDURE glob_sum_1d, glob_sum_2d, glob_sum_3d, &
& glob_sum_2d_a, glob_sum_3d_a
+ END INTERFACE
+ INTERFACE glob_sum_full
+ MODULE PROCEDURE glob_sum_full_2d, glob_sum_full_3d
END INTERFACE
INTERFACE glob_min
@@ -156,4 +160,44 @@
!
END FUNCTION glob_sum_3d_a
+
+ FUNCTION glob_sum_full_2d( ptab )
+ !!----------------------------------------------------------------------
+ !! *** FUNCTION glob_sum_2d ***
+ !!
+ !! ** Purpose : perform a sum in calling DDPDD routine
+ !!----------------------------------------------------------------------
+ REAL(wp), INTENT(in), DIMENSION(:,:) :: ptab
+ REAL(wp) :: glob_sum_full_2d ! global masked sum
+ !!
+ !!-----------------------------------------------------------------------
+ !
+ glob_sum_full_2d = SUM( ptab(:,:) )
+ IF( lk_mpp ) CALL mpp_sum( glob_sum_full_2d )
+ !
+ END FUNCTION glob_sum_full_2d
+
+ FUNCTION glob_sum_full_3d( ptab )
+ !!----------------------------------------------------------------------
+ !! *** FUNCTION glob_sum_3d ***
+ !!
+ !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine
+ !!----------------------------------------------------------------------
+ REAL(wp), INTENT(in), DIMENSION(:,:,:) :: ptab
+ REAL(wp) :: glob_sum_full_3d ! global masked sum
+ !!
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ INTEGER :: ijpk ! local variables: size of ptab
+ !!-----------------------------------------------------------------------
+ !
+ ijpk = SIZE(ptab,3)
+ !
+ glob_sum_3d = 0.e0
+ DO jk = 1, ijpk
+ glob_sum_3d = glob_sum_3d + SUM( ptab(:,:,jk) )
+ END DO
+ IF( lk_mpp ) CALL mpp_sum( glob_sum_full_3d )
+ !
+ END FUNCTION glob_sum_full_3d
+
#else
@@ -314,4 +358,65 @@
END FUNCTION glob_sum_3d_a
+ FUNCTION glob_sum_full_2d( ptab )
+ !!----------------------------------------------------------------------
+ !! *** FUNCTION glob_sum_2d ***
+ !!
+ !! ** Purpose : perform a sum in calling DDPDD routine
+ !!----------------------------------------------------------------------
+ REAL(wp), INTENT(in), DIMENSION(:,:) :: ptab
+ REAL(wp) :: glob_sum_full_2d ! global masked sum
+ !!
+ COMPLEX(wp):: ctmp
+ REAL(wp) :: ztmp
+ INTEGER :: ji, jj ! dummy loop indices
+ !!-----------------------------------------------------------------------
+ !
+ ztmp = 0.e0
+ ctmp = CMPLX( 0.e0, 0.e0, wp )
+ DO jj = 1, jpj
+ DO ji =1, jpi
+ ztmp = ptab(ji,jj)
+ CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+ END DO
+ END DO
+ IF( lk_mpp ) CALL mpp_sum( ctmp ) ! sum over the global domain
+ glob_sum_full_2d = REAL(ctmp,wp)
+ !
+ END FUNCTION glob_sum_full_2d
+
+ FUNCTION glob_sum_full_3d( ptab )
+ !!----------------------------------------------------------------------
+ !! *** FUNCTION glob_sum_3d ***
+ !!
+ !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine
+ !!----------------------------------------------------------------------
+ REAL(wp), INTENT(in), DIMENSION(:,:,:) :: ptab
+ REAL(wp) :: glob_sum_full_3d ! global masked sum
+ !!
+ COMPLEX(wp):: ctmp
+ REAL(wp) :: ztmp
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ INTEGER :: ijpk ! local variables: size of ptab
+ !!-----------------------------------------------------------------------
+ !
+ ijpk = SIZE(ptab,3)
+ !
+ ztmp = 0.e0
+ ctmp = CMPLX( 0.e0, 0.e0, wp )
+ DO jk = 1, ijpk
+ DO jj = 1, jpj
+ DO ji =1, jpi
+ ztmp = ptab(ji,jj,jk)
+ CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+ END DO
+ END DO
+ END DO
+ IF( lk_mpp ) CALL mpp_sum( ctmp ) ! sum over the global domain
+ glob_sum_full_3d = REAL(ctmp,wp)
+ !
+ END FUNCTION glob_sum_full_3d
+
+
+
#endif