Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
No results found
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 1963 additions and 263 deletions
src/OCE/ICB/icbrst.F90
View file @ 672ff863
......@@ -24,6 +24,7 @@ MODULE icbrst
USE lib_mpp ! NEMO MPI library, lk_mpp in particular
USE netcdf ! netcdf routines for IO
USE iom
USE ioipsl, ONLY : ju2ymds ! for calendar
USE icb_oce ! define iceberg arrays
USE icbutl ! iceberg utility routines
......@@ -190,9 +191,12 @@ CONTAINS
INTEGER :: jn ! dummy loop index
INTEGER :: idg ! number of digits
INTEGER :: ix_dim, iy_dim, ik_dim, in_dim
INTEGER :: iyear, imonth, iday
REAL (wp) :: zsec
REAL (wp) :: zfjulday
CHARACTER(len=256) :: cl_path
CHARACTER(len=256) :: cl_filename
CHARACTER(len=8 ) :: cl_kt
CHARACTER(LEN=20) :: cl_kt ! ocean time-step deine as a character
CHARACTER(LEN=12 ) :: clfmt ! writing format
TYPE(iceberg), POINTER :: this
TYPE(point) , POINTER :: pt
......@@ -212,8 +216,17 @@ CONTAINS
IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/'
!
! file name
WRITE(cl_kt, '(i8.8)') kt
cl_filename = TRIM(cexper)//"_"//cl_kt//"_"//TRIM(cn_icbrst_out)
IF ( ln_rstdate ) THEN
zfjulday = fjulday + rdt / rday
IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday ) zfjulday = REAL(NINT(zfjulday),wp) ! avoid truncation error
CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )
WRITE(cl_kt, '(i4.4,2i2.2)') iyear, imonth, iday
ELSE
IF( kt > 999999999 ) THEN ; WRITE(cl_kt, * ) kt
ELSE ; WRITE(cl_kt, '(i8.8)') kt
ENDIF
ENDIF
cl_filename = TRIM(cexper)//"_"//TRIM(cl_kt)//"_"//TRIM(cn_icbrst_out)
IF( lk_mpp ) THEN
idg = MAX( INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1, 4 ) ! how many digits to we need to write? min=4, max=9
WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg ! '(a,a,ix.x,a)'
......
src/OCE/IOM/in_out_manager.F90
View file @ 672ff863
......@@ -27,6 +27,7 @@ MODULE in_out_manager
CHARACTER(lc) :: cn_ocerst_outdir !: restart output directory
LOGICAL :: ln_rstart !: start from (F) rest or (T) a restart file
LOGICAL :: ln_rst_list !: output restarts at list of times (T) or by frequency (F)
LOGICAL :: ln_rst_eos !: check equation of state used for the restart is consistent with model
INTEGER :: nn_rstctl !: control of the time step (0, 1 or 2)
INTEGER :: nn_rstssh = 0 !: hand made initilization of ssh or not (1/0)
INTEGER :: nn_it000 !: index of the first time step
......@@ -39,6 +40,7 @@ MODULE in_out_manager
INTEGER :: nn_stock !: restart file frequency
INTEGER, DIMENSION(10) :: nn_stocklist !: restart dump times
LOGICAL :: ln_mskland !: mask land points in NetCDF outputs (costly: + ~15%)
LOGICAL :: ln_rstdate !: T=> stamp output restart files with date instead of timestep
LOGICAL :: ln_cfmeta !: output additional data to netCDF files required for compliance with the CF metadata standard
LOGICAL :: ln_clobber !: clobber (overwrite) an existing file
INTEGER :: nn_chunksz !: chunksize (bytes) for NetCDF file (works only with iom_nf90 routines)
......
src/OCE/IOM/restart.F90
View file @ 672ff863
......@@ -39,6 +39,7 @@ MODULE restart
!
USE in_out_manager ! I/O manager
USE iom ! I/O module
USE ioipsl, ONLY : ju2ymds ! for calendar
USE lib_mpp ! distribued memory computing library
IMPLICIT NONE
......@@ -72,6 +73,9 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time-step
!!
INTEGER :: iyear, imonth, iday
REAL (wp) :: zsec
REAL (wp) :: zfjulday
CHARACTER(LEN=20) :: clkt ! ocean time-step deine as a character
CHARACTER(LEN=50) :: clname ! ocean output restart file name
CHARACTER(lc) :: clpath ! full path to ocean output restart file
......@@ -103,8 +107,16 @@ CONTAINS
IF( kt == nitrst - 1 .OR. nn_stock == 1 .OR. ( kt == nitend .AND. .NOT. lrst_oce ) ) THEN
IF( nitrst <= nitend .AND. nitrst > 0 ) THEN
! beware of the format used to write kt (default is i8.8, that should be large enough...)
IF( nitrst > 999999999 ) THEN ; WRITE(clkt, * ) nitrst
ELSE ; WRITE(clkt, '(i8.8)') nitrst
IF ( ln_rstdate ) THEN
zfjulday = fjulday + rdt / rday
IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday ) zfjulday = REAL(NINT(zfjulday),wp) ! avoid truncation error
CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )
WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
ELSE
IF( nitrst > 999999999 ) THEN ; WRITE(clkt, * ) nitrst
ELSE ; WRITE(clkt, '(i8.8)') nitrst
ENDIF
ENDIF
! create the file
clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_ocerst_out)
......@@ -183,6 +195,9 @@ CONTAINS
#endif
ENDIF
CALL iom_rstput( kt, nitrst, numrow, 'neos' , REAL(neos)) ! equation of state
!CALL iom_rstput( kt, nitrst, numrow, 'neos' , neos , ktype = jp_i1, ldxios = lwxios) ! equation of state
IF( ln_diurnal ) CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst )
IF( kt == nitrst ) THEN
IF( .NOT.lwxios ) THEN
......@@ -267,6 +282,8 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices
INTEGER :: jk
INTEGER :: id1
REAL(wp) :: zeos
REAL(wp), DIMENSION(jpi, jpj, jpk) :: w3d
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zgdept ! 3D workspace for QCO
!!----------------------------------------------------------------------
......@@ -299,6 +316,18 @@ CONTAINS
CALL iom_get( numror, jpdom_auto, 'tn', ts(:,:,:,jp_tem,Kmm) )
CALL iom_get( numror, jpdom_auto, 'sn', ts(:,:,:,jp_sal,Kmm) )
!
IF ( ln_rst_eos ) THEN
! Check equation of state used is consistent with the restart
IF( iom_varid( numror, 'neos') == -1 ) THEN
CALL ctl_stop( 'restart, rst_read: variable neos not found. STOP check that the equations of state in the restart file and in the namelist nameos are consistent and use ln_rst_eos=F')
ELSE
CALL iom_get( numror, 'neos' , zeos )
IF ( INT(zeos) /= neos ) CALL ctl_stop( 'restart, rst_read: equation of state used in restart file differs from namelist nameos')
ENDIF
ENDIF
IF( l_1st_euler ) THEN !* Euler restart (MLF only)
IF(lwp) WRITE(numout,*) ' Kbb u, v and T-S fields set to Kmm values'
uu(:,:,: ,Kbb) = uu(:,:,: ,Kmm) ! all before fields set to now values
......
src/OCE/LBC/mppini.F90
View file @ 672ff863
......@@ -333,8 +333,8 @@ CONTAINS
nimpp = iimppt(ii,ij)
njmpp = ijmppt(ii,ij)
!
CALL init_doloop ! set start/end indices of do-loop, depending on the halo width value (nn_hls)
CALL init_locglo ! define now functions needed to convert indices from/to global to/from local domains
CALL init_doloop ! set start/end indices of do-loop, depending on the halo width value (nn_hls)
!
IF(lwp) THEN
WRITE(numout,*)
......@@ -1409,9 +1409,46 @@ ENDIF
Njs0 = 1+nn_hls
Nie0 = jpi-nn_hls
Nje0 = jpj-nn_hls
!
Ni_0 = Nie0 - Nis0 + 1
Nj_0 = Nje0 - Njs0 + 1
! Start and end indexes for actual coupling fields on extended grid
Nis0_ext = Nis0
Nie0_ext = Nie0
Njs0_ext = Njs0
Nje0_ext = Nje0
IF (mig(1) == 1 ) THEN
! Drag start column 1 place to the left/west
Nis0_ext=Nis0_ext-1
ENDIF
IF (mig(jpi) == jpiglo ) THEN
! Drag end column 1 place to the right/east
Nie0_ext=Nie0_ext+1
ENDIF
! RSRH we don't adjust anything in the N-S (j) direction
! since the content of the N-fold is catared for by populating values
! using lbc_lnk rather than relying on getting them from the coupler.
! This is rather confused by the fact that jpjglo has a value BIGGER
! than it did at pre 4.2... e.g. for ORCA1 it's set to 333 instead of 332
! which is rather baffling and confuses some dimensioning calculations
! which previously worked fine pre 4.2.
! Set up dimensions for old style coupling exchanges on extended grid
Ni_0_ext = Ni_0
Nj_0_ext = Nj_0
IF (mig(1) == 1 .OR. mig(jpi) == jpiglo ) THEN
! We're at the extreme left or right edge of the grid so need to cater
! for an extra column
Ni_0_ext = Ni_0_ext + 1
ENDIF
!
jpkm1 = jpk-1 ! " "
!
......
src/OCE/LDF/ldftra.F90
View file @ 672ff863
......@@ -71,6 +71,8 @@ MODULE ldftra
INTEGER , PUBLIC :: nn_aei_ijk_t !: choice of time/space variation of the eiv coeff.
REAL(wp), PUBLIC :: rn_Ue !: lateral diffusive velocity [m/s]
REAL(wp), PUBLIC :: rn_Le !: lateral diffusive length [m]
INTEGER, PUBLIC :: nn_ldfeiv_shape !: shape of bounding coefficient (Treguier et al formulation only)
! ! Flag to control the type of lateral diffusive operator
INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 ! error in specification of lateral diffusion
......@@ -495,7 +497,8 @@ CONTAINS
REAL(wp) :: zah_max, zUfac ! - scalar
!!
NAMELIST/namtra_eiv/ ln_ldfeiv , ln_ldfeiv_dia, & ! eddy induced velocity (eiv)
& nn_aei_ijk_t, rn_Ue, rn_Le ! eiv coefficient
& nn_aei_ijk_t, rn_Ue, rn_Le, & ! eiv coefficient
& nn_ldfeiv_shape
!!----------------------------------------------------------------------
!
IF(lwp) THEN ! control print
......@@ -588,6 +591,7 @@ CONTAINS
IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( latitude, longitude, time )'
IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )'
IF(lwp) WRITE(numout,*) ' maximum allowed value: aei0 = ', aei0, ' m2/s'
IF(lwp) WRITE(numout,*) ' shape of bounding coefficient : ',nn_ldfeiv_shape
!
l_ldfeiv_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90
!
......@@ -636,14 +640,21 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: paeiu, paeiv ! eiv coefficient [m2/s]
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zfw, ze3w, zn2, z1_f20, zzaei ! local scalars
REAL(wp), DIMENSION(jpi,jpj) :: zn, zah, zhw, zRo, zaeiw ! 2D workspace
REAL(wp) :: zfw, ze3w, zn2, z1_f20, zzaei, z2_3 ! local scalars
REAL(wp), DIMENSION(jpi,jpj) :: zn, zah, zhw, zRo, zRo_lim, zTclinic_recip, zaeiw, zratio ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmodslp ! 3D workspace
!!----------------------------------------------------------------------
!
zn (:,:) = 0._wp ! Local initialization
zmodslp(:,:,:) = 0._wp
zhw(:,:) = 5._wp
zah(:,:) = 0._wp
zRo(:,:) = 0._wp
zRo_lim(:,:) = 0._wp
zTclinic_recip(:,:) = 0._wp
zratio(:,:) = 0._wp
zaeiw(:,:) = 0._wp
! ! Compute lateral diffusive coefficient at T-point
IF( ln_traldf_triad ) THEN
DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
......@@ -670,8 +681,9 @@ CONTAINS
! eddies using the isopycnal slopes calculated in ldfslp.F :
! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
& + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w
zmodslp(ji,jj,jk) = wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
& + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)
zah(ji,jj) = zah(ji,jj) + zn2 * zmodslp(ji,jj,jk) * ze3w
zhw(ji,jj) = zhw(ji,jj) + ze3w
END_3D
ENDIF
......@@ -679,17 +691,69 @@ CONTAINS
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 )
! Rossby radius at w-point taken betwenn 2 km and 40km
zRo(ji,jj) = MAX( 2.e3 , MIN( .4 * zn(ji,jj) / zfw, 40.e3 ) )
zRo(ji,jj) = .4 * zn(ji,jj) / zfw
zRo_lim(ji,jj) = MAX( 2.e3 , MIN( zRo(ji,jj), 40.e3 ) )
! Compute aeiw by multiplying Ro^2 and T^-1
zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)
zTclinic_recip(ji,jj) = SQRT( MAX(zah(ji,jj),0._wp) / zhw(ji,jj) ) * tmask(ji,jj,1)
zaeiw(ji,jj) = zRo_lim(ji,jj) * zRo_lim(ji,jj) * zTclinic_recip(ji,jj)
END_2D
IF( iom_use('N_2d') ) CALL iom_put('N_2d',zn(:,:)/zhw(:,:))
IF( iom_use('modslp') ) CALL iom_put('modslp',SQRT(zmodslp(:,:,:)) )
CALL iom_put('RossRad',zRo)
CALL iom_put('RossRadlim',zRo_lim)
CALL iom_put('Tclinic_recip',zTclinic_recip)
! !== Bound on eiv coeff. ==!
z1_f20 = 1._wp / ( 2._wp * omega * sin( rad * 20._wp ) )
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease
zaeiw(ji,jj) = MIN( zzaei , paei0 ) ! Max value = paei0
END_2D
z2_3 = 2._wp/3._wp
SELECT CASE(nn_ldfeiv_shape)
CASE(0) !! Standard shape applied - decrease in tropics and cap.
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease
zaeiw(ji,jj) = MIN( zzaei, paei0 )
END_2D
CASE(1) !! Abrupt cut-off on Rossby radius:
! JD : modifications here to introduce scaling by local rossby radius of deformation vs local grid scale
! arbitrary decision that GM is de-activated if local rossy radius larger than 2 times local grid scale
! based on Hallberg (2013)
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
IF ( zRo(ji,jj) >= ( 2._wp * MIN( e1t(ji,jj), e2t(ji,jj) ) ) ) THEN
! TODO : use a version of zRo that integrates over a few time steps ?
zaeiw(ji,jj) = 0._wp
ELSE
zaeiw(ji,jj) = MIN( zaeiw(ji,jj), paei0 )
ENDIF
END_2D
CASE(2) !! Rossby radius ramp type 1:
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zratio(ji,jj) = zRo(ji,jj)/MIN(e1t(ji,jj),e2t(ji,jj))
zaeiw(ji,jj) = MIN( zaeiw(ji,jj), MAX( 0._wp, MIN( 1._wp, z2_3*(2._wp - zratio(ji,jj)) ) ) * paei0 )
END_2D
CALL iom_put('RR_GS',zratio)
CASE(3) !! Rossby radius ramp type 2:
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zratio(ji,jj) = MIN(e1t(ji,jj),e2t(ji,jj))/zRo(ji,jj)
zaeiw(ji,jj) = MIN( zaeiw(ji,jj), MAX( 0._wp, MIN( 1._wp, z2_3*( zratio(ji,jj) - 0.5_wp ) ) ) * paei0 )
END_2D
CASE(4) !! bathymetry ramp:
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zaeiw(ji,jj) = MIN( zaeiw(ji,jj), MAX( 0._wp, MIN( 1._wp, 0.001*(ht_0(ji,jj) - 2000._wp) ) ) * paei0 )
END_2D
CASE(5) !! Rossby radius ramp type 1 applied to Treguier et al coefficient rather than cap:
!! Note the ramp is RR/GS=[2.0,1.0] (not [2.0,0.5] as for cases 2,3) and we ramp up
!! to 5% of the Treguier et al coefficient, aiming for peak values of around 100m2/s
!! at high latitudes rather than 2000m2/s which is what you get in eORCA025 with an
!! uncapped coefficient.
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zratio(ji,jj) = zRo(ji,jj)/MIN(e1t(ji,jj),e2t(ji,jj))
zaeiw(ji,jj) = MAX( 0._wp, MIN( 1._wp, 2._wp - zratio(ji,jj) ) ) * 0.05 * zaeiw(ji,jj)
zaeiw(ji,jj) = MIN( zaeiw(ji,jj), paei0 )
END_2D
CALL iom_put('RR_GS',zratio)
CASE DEFAULT
CALL ctl_stop('ldf_eiv: Unrecognised option for nn_ldfeiv_shape.')
END SELECT
IF( nn_hls == 1 ) CALL lbc_lnk( 'ldftra', zaeiw(:,:), 'W', 1.0_wp ) ! lateral boundary condition
!
DO_2D( 0, 0, 0, 0 )
......
src/OCE/SBC/cpl_oasis3.F90
View file @ 672ff863
......@@ -30,20 +30,35 @@ MODULE cpl_oasis3
USE xios ! I/O server
#endif
USE par_oce ! ocean parameters
USE cpl_rnf_1d, ONLY: nn_cpl_river ! Variables used in 1D river outflow
USE dom_oce ! ocean space and time domain
USE in_out_manager ! I/O manager
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE lib_mpp
#if defined key_agrif || ! defined key_mpi_off
USE MPI
#endif
IMPLICIT NONE
PRIVATE
#if ! defined key_oasis3
! Dummy interface to oasis_get if not using oasis
INTERFACE oasis_get
MODULE PROCEDURE oasis_get_1d, oasis_get_2d
END INTERFACE
#endif
PUBLIC cpl_init
PUBLIC cpl_define
PUBLIC cpl_snd
PUBLIC cpl_rcv
PUBLIC cpl_rcv_1d
PUBLIC cpl_freq
PUBLIC cpl_finalize
INTEGER, PARAMETER :: localRoot = 0
INTEGER, PUBLIC :: OASIS_Rcv = 1 !: return code if received field
INTEGER, PUBLIC :: OASIS_idle = 0 !: return code if nothing done by oasis
INTEGER :: ncomp_id ! id returned by oasis_init_comp
......@@ -67,9 +82,13 @@ MODULE cpl_oasis3
INTEGER :: nrcv ! total number of fields received
INTEGER :: nsnd ! total number of fields sent
INTEGER :: ncplmodel ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
INTEGER, PUBLIC, PARAMETER :: nmaxfld=62 ! Maximum number of coupling fields
INTEGER, PUBLIC, PARAMETER :: nmaxfld=100 ! Maximum number of coupling fields
INTEGER, PUBLIC, PARAMETER :: nmaxcat=5 ! Maximum number of coupling fields
INTEGER, PUBLIC, PARAMETER :: nmaxcpl=5 ! Maximum number of coupling fields
LOGICAL, PARAMETER :: ltmp_wapatch = .FALSE. ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define
LOGICAL, PARAMETER :: ltmp_landproc = .TRUE. ! patch to restrict coupled area to non halo cells
TYPE, PUBLIC :: FLD_CPL !: Type for coupling field information
LOGICAL :: laction ! To be coupled or not
......@@ -79,11 +98,16 @@ MODULE cpl_oasis3
INTEGER, DIMENSION(nmaxcat,nmaxcpl) :: nid ! Id of the field (no more than 9 categories and 9 extrena models)
INTEGER :: nct ! Number of categories in field
INTEGER :: ncplmodel ! Maximum number of models to/from which this variable may be sent/received
INTEGER :: dimensions ! Number of dimensions of coupling field
END TYPE FLD_CPL
TYPE(FLD_CPL), DIMENSION(nmaxfld), PUBLIC :: srcv, ssnd !: Coupling fields
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: exfld ! Temporary buffer for receiving
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: exfld_ext ! Temporary buffer for receiving with wrap points
INTEGER :: ishape_ext(4) ! shape of 2D arrays passed to PSMILe extended for wrap points in weights data
INTEGER :: ishape(4) ! shape of 2D arrays passed to PSMILe
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
......@@ -103,6 +127,7 @@ CONTAINS
!!--------------------------------------------------------------------
CHARACTER(len = *), INTENT(in ) :: cd_modname ! model name as set in namcouple file
INTEGER , INTENT( out) :: kl_comm ! local communicator of the model
INTEGER :: error
!!--------------------------------------------------------------------
! WARNING: No write in numout in this routine
......@@ -123,6 +148,7 @@ CONTAINS
IF( nerror /= OASIS_Ok ) &
CALL oasis_abort (ncomp_id, 'cpl_init','Failure in oasis_get_localcomm' )
!
END SUBROUTINE cpl_init
......@@ -138,13 +164,26 @@ CONTAINS
INTEGER, INTENT(in) :: krcv, ksnd ! Number of received and sent coupling fields
INTEGER, INTENT(in) :: kcplmodel ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
!
INTEGER :: id_part
INTEGER :: id_part_0d ! Partition for 0d fields
INTEGER :: id_part_rnf_1d ! Partition for 1d river outflow fields
INTEGER :: id_part_2d ! Partition for 2d fields
INTEGER :: id_part_2d_ext ! Partition for 2d fields extended for old (pre vn4.2) style remapping weights!
INTEGER :: id_part_temp ! Temperary partition used to choose either 0d or 1d partitions
INTEGER :: paral(5) ! OASIS3 box partition
INTEGER :: ishape(4) ! shape of arrays passed to PSMILe
INTEGER :: paral_ext(5) ! OASIS3 box partition extended
INTEGER :: ishape0d1d(2) ! Shape of 0D or 1D arrays passed to PSMILe.
INTEGER :: var_nodims(2) ! Number of coupling field dimensions.
! var_nodims(1) is redundant from OASIS3-MCT vn4.0 onwards
! but retained for backward compatibility.
! var_nodims(2) is the number of fields in a bundle
! or 1 for unbundled fields (bundles are not yet catered for
! in NEMO hence we default to 1).
INTEGER :: ji,jc,jm ! local loop indicees
LOGICAL :: llenddef ! should we call xios_oasis_enddef and oasis_enddef?
CHARACTER(LEN=64) :: zclname
CHARACTER(LEN=2) :: cli2
INTEGER :: i_offset ! Used in calculating offset for extended partition.
!!--------------------------------------------------------------------
IF(lwp) WRITE(numout,*)
......@@ -173,10 +212,13 @@ CONTAINS
ishape(2) = Ni_0
ishape(3) = 1
ishape(4) = Nj_0
!
ishape0d1d(1) = 0
ishape0d1d(2) = 0 !
! ... Allocate memory for data exchange
!
ALLOCATE(exfld(Ni_0, Nj_0), stat = nerror) ! allocate only inner domain (without halos)
ALLOCATE(exfld(Ni_0, Nj_0), stat = nerror) ! allocate full domain (without halos)
IF( nerror > 0 ) THEN
CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in allocating exfld') ; RETURN
ENDIF
......@@ -198,7 +240,96 @@ CONTAINS
WRITE(numout,*) ' multiexchg: Njs0, Nje0, njmpp =', Njs0, Nje0, njmpp
ENDIF
CALL oasis_def_partition ( id_part, paral, nerror, Ni0glo*Nj0glo ) ! global number of points, excluding halos
! We still set up the new vn4.2 style box partition for reference, though it doesn't actually get used,
! we can easily swap back to it if we ever manage to successfully generate vn4.2 compatible weights, or introduce
! RTL controls to distinguish between onl and new style weights.
CALL oasis_def_partition ( id_part_2d, paral, nerror, Ni0glo*Nj0glo ) ! global number of points, excluding halos
! RSRH Set up 2D box partition for compatibility with existing weights files
! so we don't have to generate and manage multiple sets of weights purely because of
! the changes to nemo 4.2+ code!
! This is just a hack for global cyclic models for the time being
Ni0glo_ext = jpiglo
Nj0glo_ext = Nj0glo +1 ! We can't use jpjglo here because for some reason at 4.2 this is bigger
! than at 4.0.... e.g. for ORCA1 it is 333 when it should only be 332!
! RSRH extended shapes for old style dimensioning. Allows backwards compatibility with existing weights files,
! which the new code DOES NOT, causing headaches not only for users but also for management of weights files.
ishape_ext(:) = ishape(:)
IF (mig(1) == 1 .OR. mig(jpi)==jpiglo) THEN
! Extra columns in PEs dealing with wrap points
ishape_ext(2) = ishape_ext(2) + 1
ENDIF
! Workout any extra offset in the i dimension
IF (mig(1) == 1 ) THEN
i_offset = mig0(nn_hls)
ELSE
i_offset = mig(nn_hls)
ENDIF
ALLOCATE(exfld_ext(ishape_ext(2), ishape_ext(4)), stat = nerror) ! allocate full domain (with wrap pts)
IF( nerror > 0 ) THEN
CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in allocating exfld_ext') ; RETURN
ENDIF
! Now we have the appropriate dimensions, we can set up the partition array for the old-style extended grid
paral_ext(1) = 2 ! box partitioning
paral_ext(2) = (Ni0glo_ext * mjg0(nn_hls)) + i_offset ! NEMO lower left corner global offset, with wrap pts
paral_ext(3) = Ni_0_ext ! local extent in i, including halos
paral_ext(4) = Nj_0_ext ! local extent in j, including halos
paral_ext(5) = Ni0glo_ext ! global extent in x, including halos
IF( sn_cfctl%l_oasout ) THEN
WRITE(numout,*) ' multiexchg: paral_ext (1:5)', paral_ext, jpiglo, jpjglo, Ni0glo_ext, Nj0glo_ext
WRITE(numout,*) ' multiexchg: Ni_0_ext, Nj_0_ext i_offset =', Ni_0_ext, Nj_0_ext, i_offset
WRITE(numout,*) ' multiexchg: Nis0_ext, Nie0_ext =', Nis0_ext, Nie0_ext
WRITE(numout,*) ' multiexchg: Njs0_ext, Nje0_ext =', Njs0_ext, Nje0_ext
ENDIF
! Define our extended grid
CALL oasis_def_partition ( id_part_2d_ext, paral_ext, nerror, Ni0glo_ext*Nj0glo_ext )
! OK so now we should have a usable 2d partition for fields defined WITH redundant points.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! A special partition is needed for 0D fields
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
paral(1) = 0 ! serial partitioning
paral(2) = 0
IF ( mpprank == 0) THEN
paral(3) = 1 ! Size of array to couple (scalar)
ELSE
paral(3) = 0 ! Dummy size for PE's not involved
END IF
paral(4) = 0
paral(5) = 0
CALL oasis_def_partition ( id_part_0d, paral, nerror, 1 )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Another special partition is needed for 1D river routing fields
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
paral(1) = 0 ! serial partitioning
paral(2) = 0
IF ( mpprank == 0) THEN
paral(3) = nn_cpl_river ! Size of array to couple (vector)
ELSE
paral(3) = 0 ! Dummy size for PE's not involved
END IF
paral(4) = 0
paral(5) = 0
CALL oasis_def_partition ( id_part_rnf_1d, paral, nerror, nn_cpl_river )
!
! ... Announce send variables.
!
......@@ -232,14 +363,24 @@ CONTAINS
ENDIF
#endif
IF( sn_cfctl%l_oasout ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_Out
CALL oasis_def_var (ssnd(ji)%nid(jc,jm), zclname, id_part , (/ 2, 1 /), &
& OASIS_Out , ishape , OASIS_REAL, nerror )
flush(numout)
IF( sn_cfctl%l_oasout ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_Out
CALL oasis_def_var (ssnd(ji)%nid(jc,jm), zclname, id_part_2d_ext , (/ 2, 1 /), &
& OASIS_Out , ishape_ext , OASIS_REAL, nerror )
IF( nerror /= OASIS_Ok ) THEN
WRITE(numout,*) 'Failed to define transient ', ji, jc, jm, " "//TRIM(zclname)
CALL oasis_abort ( ssnd(ji)%nid(jc,jm), 'cpl_define', 'Failure in oasis_def_var' )
ENDIF
IF( sn_cfctl%l_oasout .AND. ssnd(ji)%nid(jc,jm) /= -1 ) WRITE(numout,*) "variable defined in the namcouple"
IF( sn_cfctl%l_oasout .AND. ssnd(ji)%nid(jc,jm) == -1 ) WRITE(numout,*) "variable NOT defined in the namcouple"
IF( sn_cfctl%l_oasout .AND. ssnd(ji)%nid(jc,jm) /= -1 ) WRITE(numout,*) "put variable defined in the namcouple"
IF( sn_cfctl%l_oasout .AND. ssnd(ji)%nid(jc,jm) == -1 ) WRITE(numout,*) "put variable NOT defined in the namcouple"
END DO
END DO
ENDIF
......@@ -276,15 +417,57 @@ CONTAINS
zclname=TRIM(Agrif_CFixed())//'_'//TRIM(zclname)
ENDIF
#endif
IF( sn_cfctl%l_oasout ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_In
CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part , (/ 2, 1 /), &
& OASIS_In , ishape , OASIS_REAL, nerror )
flush(numout)
! Define 0D (Greenland or Antarctic ice mass) or 1D (river outflow) coupling fields
IF (srcv(ji)%dimensions <= 1) THEN
var_nodims(1) = 1
var_nodims(2) = 1 ! Modify this value to cater for bundled fields.
IF (mpprank == 0) THEN
IF (srcv(ji)%dimensions == 0) THEN
! If 0D then set temporary variables to 0D components
id_part_temp = id_part_0d
ishape0d1d(2) = 1
ELSE
! If 1D then set temporary variables to river outflow components
id_part_temp = id_part_rnf_1d
ishape0d1d(2)= nn_cpl_river
END IF
CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part_temp , var_nodims, &
OASIS_In , ishape0d1d(1:2) , OASIS_REAL, nerror )
ELSE
! Dummy call to keep OASIS3-MCT happy.
CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part_0d , var_nodims, &
OASIS_In , ishape0d1d(1:2) , OASIS_REAL, nerror )
END IF
ELSE
! It's a "normal" 2D (or pseudo 3D) coupling field.
! ... Set the field dimension and bundle count
var_nodims(1) = 2
var_nodims(2) = 1 ! Modify this value to cater for bundled fields.
CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part_2d_ext , var_nodims, &
OASIS_In , ishape_ext , OASIS_REAL, nerror )
ENDIF
IF( nerror /= OASIS_Ok ) THEN
WRITE(numout,*) 'Failed to define transient ', ji, jc, jm, " "//TRIM(zclname)
CALL oasis_abort ( srcv(ji)%nid(jc,jm), 'cpl_define', 'Failure in oasis_def_var' )
ENDIF
IF( sn_cfctl%l_oasout .AND. srcv(ji)%nid(jc,jm) /= -1 ) WRITE(numout,*) "variable defined in the namcouple"
IF( sn_cfctl%l_oasout .AND. srcv(ji)%nid(jc,jm) == -1 ) WRITE(numout,*) "variable NOT defined in the namcouple"
IF( sn_cfctl%l_oasout .AND. srcv(ji)%nid(jc,jm) /= -1 ) WRITE(numout,*) "get variable defined in the namcouple"
IF( sn_cfctl%l_oasout .AND. srcv(ji)%nid(jc,jm) == -1 ) WRITE(numout,*) "get variable NOT defined in the namcouple"
END DO
END DO
......@@ -330,13 +513,17 @@ CONTAINS
INTEGER :: jc,jm ! local loop index
!!--------------------------------------------------------------------
!
! snd data to OASIS3
!
DO jc = 1, ssnd(kid)%nct
DO jm = 1, ssnd(kid)%ncplmodel
IF( ssnd(kid)%nid(jc,jm) /= -1 ) THEN ! exclude halos from data sent to oasis
CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(Nis0:Nie0, Njs0:Nje0,jc), kinfo )
! The field is "put" directly, using appropriate start/end indexing - i.e. we don't
! copy it to an intermediate buffer.
CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc), kinfo )
IF ( sn_cfctl%l_oasout ) THEN
IF ( kinfo == OASIS_Sent .OR. kinfo == OASIS_ToRest .OR. &
......@@ -346,10 +533,11 @@ CONTAINS
WRITE(numout,*) 'oasis_put: ivarid ', ssnd(kid)%nid(jc,jm)
WRITE(numout,*) 'oasis_put: kstep ', kstep
WRITE(numout,*) 'oasis_put: info ', kinfo
WRITE(numout,*) ' - Minimum value is ', MINVAL(pdata(Nis0:Nie0,Njs0:Nje0,jc))
WRITE(numout,*) ' - Maximum value is ', MAXVAL(pdata(Nis0:Nie0,Njs0:Nje0,jc))
WRITE(numout,*) ' - Sum value is ', SUM(pdata(Nis0:Nie0,Njs0:Nje0,jc))
WRITE(numout,*) ' - Minimum value is ', MINVAL(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
WRITE(numout,*) ' - Maximum value is ', MAXVAL(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
WRITE(numout,*) ' - Sum value is ', SUM(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
WRITE(numout,*) '****************'
CALL FLUSH(numout)
ENDIF
ENDIF
......@@ -357,6 +545,7 @@ CONTAINS
ENDDO
ENDDO
!
END SUBROUTINE cpl_snd
......@@ -375,13 +564,16 @@ CONTAINS
INTEGER , INTENT( out) :: kinfo ! OASIS3 info argument
!!
INTEGER :: jc,jm ! local loop index
LOGICAL :: llaction, ll_1st
LOGICAL :: llaction, ll_1st, lrcv
!!--------------------------------------------------------------------
!
! receive local data from OASIS3 on every process
!
kinfo = OASIS_idle
lrcv=.FALSE.
!
DO jc = 1, srcv(kid)%nct
ll_1st = .TRUE.
......@@ -389,7 +581,7 @@ CONTAINS
IF( srcv(kid)%nid(jc,jm) /= -1 ) THEN
CALL oasis_get ( srcv(kid)%nid(jc,jm), kstep, exfld, kinfo )
CALL oasis_get ( srcv(kid)%nid(jc,jm), kstep, exfld_ext, kinfo )
llaction = kinfo == OASIS_Recvd .OR. kinfo == OASIS_FromRest .OR. &
& kinfo == OASIS_RecvOut .OR. kinfo == OASIS_FromRestOut
......@@ -398,14 +590,18 @@ CONTAINS
& WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , llaction, kinfo, kstep, srcv(kid)%nid(jc,jm)
IF( llaction ) THEN ! data received from oasis do not include halos
! but DO still cater for wrap columns when using pre vn4.2 compatible remapping weights.
kinfo = OASIS_Rcv
lrcv=.TRUE.
IF( ll_1st ) THEN
pdata(Nis0:Nie0,Njs0:Nje0,jc) = exfld(:,:) * pmask(Nis0:Nie0,Njs0:Nje0,jm)
pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc) = exfld_ext(:,:) * pmask(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jm)
ll_1st = .FALSE.
ELSE
pdata(Nis0:Nie0,Njs0:Nje0,jc) = pdata(Nis0:Nie0,Njs0:Nje0,jc) &
& + exfld(:,:) * pmask(Nis0:Nie0,Njs0:Nje0,jm)
pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc) = pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc) &
& + exfld_ext(:,:) * pmask(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jm)
ENDIF
IF ( sn_cfctl%l_oasout ) THEN
......@@ -414,10 +610,11 @@ CONTAINS
WRITE(numout,*) 'oasis_get: ivarid ' , srcv(kid)%nid(jc,jm)
WRITE(numout,*) 'oasis_get: kstep', kstep
WRITE(numout,*) 'oasis_get: info ', kinfo
WRITE(numout,*) ' - Minimum value is ', MINVAL(pdata(Nis0:Nie0,Njs0:Nje0,jc))
WRITE(numout,*) ' - Maximum value is ', MAXVAL(pdata(Nis0:Nie0,Njs0:Nje0,jc))
WRITE(numout,*) ' - Sum value is ', SUM(pdata(Nis0:Nie0,Njs0:Nje0,jc))
WRITE(numout,*) ' - Minimum value is ', MINVAL(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
WRITE(numout,*) ' - Maximum value is ', MAXVAL(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
WRITE(numout,*) ' - Sum value is ', SUM(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
WRITE(numout,*) '****************'
CALL FLUSH(numout)
ENDIF
ENDIF
......@@ -426,15 +623,130 @@ CONTAINS
ENDDO
!--- we must call lbc_lnk to fill the halos that where not received.
IF( .NOT. ll_1st ) THEN
CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,jc), srcv(kid)%clgrid, srcv(kid)%nsgn )
ENDIF
ENDDO
! RSRH I've changed this since:
! 1) it seems multi cat fields may not be properly updated in halos when called on a per
! category basis(?)
! 2) it's more efficient to have a single call (and simpler to understand) to update ALL
! categories at the same time!
!--- Call lbc_lnk to populate halos of received fields.
IF (lrcv) then
CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,:), srcv(kid)%clgrid, srcv(kid)%nsgn )
endif
!
END SUBROUTINE cpl_rcv
SUBROUTINE cpl_rcv_1d( kid, kstep, pdata, nitems, kinfo )
!!---------------------------------------------------------------------
!! *** ROUTINE cpl_rcv_1d ***
!!
!! ** Purpose : - A special version of cpl_rcv to deal exclusively with
!! receipt of 0D or 1D fields.
!! The fields are recieved into a 1D array buffer which is simply a
!! dynamically sized sized array (which may be of size 1)
!! of 0 dimensional fields. This allows us to pass miltiple 0D
!! fields via a single put/get operation.
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: nitems ! Number of 0D items to recieve
! during this get operation. i.e.
! The size of the 1D array in which
! 0D items are passed.
INTEGER , INTENT(in ) :: kid ! ID index of the incoming
! data.
INTEGER , INTENT(in ) :: kstep ! ocean time-step in seconds
REAL(wp), INTENT(inout) :: pdata(1:nitems) ! The original value(s),
! unchanged if nothing is
! received
INTEGER , INTENT( out) :: kinfo ! OASIS3 info argument
!!
REAL(wp) :: recvfld(1:nitems) ! Local receive field buffer
INTEGER :: jc,jm ! local loop index
INTEGER :: ierr
LOGICAL :: llaction
INTEGER :: MPI_WORKING_PRECISION
INTEGER :: number_to_print
!!--------------------------------------------------------------------
!
! receive local data from OASIS3 on every process
!
kinfo = OASIS_idle
!
! 0D and 1D fields won't have categories or any other form of "pseudo level"
! so we only cater for a single set of values and thus don't bother
! with a loop over the jc index
jc = 1
DO jm = 1, srcv(kid)%ncplmodel
IF( srcv(kid)%nid(jc,jm) /= -1 ) THEN
IF ( ( srcv(kid)%dimensions <= 1) .AND. (mpprank == 0) ) THEN
! Since there is no concept of data decomposition for zero
! dimension fields, they must only be exchanged through the master PE,
! unlike "normal" 2D field cases where every PE is involved.
CALL oasis_get ( srcv(kid)%nid(jc,jm), kstep, recvfld, kinfo )
llaction = kinfo == OASIS_Recvd .OR. kinfo == OASIS_FromRest .OR. &
kinfo == OASIS_RecvOut .OR. kinfo == OASIS_FromRestOut
IF ( sn_cfctl%l_oasout ) WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , &
llaction, kinfo, kstep, srcv(kid)%nid(jc,jm)
IF ( llaction ) THEN
kinfo = OASIS_Rcv
pdata(1:nitems) = recvfld(1:nitems)
IF( sn_cfctl%l_oasout ) THEN
number_to_print = 10
IF ( nitems < number_to_print ) number_to_print = nitems
WRITE(numout,*) '****************'
WRITE(numout,*) 'oasis_get: Incoming ', srcv(kid)%clname
WRITE(numout,*) 'oasis_get: ivarid ' , srcv(kid)%nid(jc,jm)
WRITE(numout,*) 'oasis_get: kstep', kstep
WRITE(numout,*) 'oasis_get: info ', kinfo
WRITE(numout,*) ' - Minimum Value is ', MINVAL(pdata(:))
WRITE(numout,*) ' - Maximum value is ', MAXVAL(pdata(:))
WRITE(numout,*) ' - Start of data is ', pdata(1:number_to_print)
WRITE(numout,*) '****************'
CALL FLUSH(numout)
ENDIF
ENDIF
ENDIF
ENDIF
ENDDO
#if defined key_mpi_off
CALL oasis_abort( ncomp_id, "cpl_rcv_1d", "Unable to use mpi_bcast with key_mpi_off in your list of NEMO keys." )
#else
! Set the precision that we want to broadcast using MPI_BCAST
SELECT CASE( wp )
CASE( sp )
MPI_WORKING_PRECISION = MPI_REAL ! Single precision
CASE( dp )
MPI_WORKING_PRECISION = MPI_DOUBLE_PRECISION ! Double precision
CASE default
CALL oasis_abort( ncomp_id, "cpl_rcv_1d", "Could not find precision for coupling 0d or 1d field" )
END SELECT
! We have to broadcast (potentially) received values from PE 0 to all
! the others. If no new data has been received we're just
! broadcasting the existing values but there's no more efficient way
! to deal with that w/o NEMO adopting a UM-style test mechanism
! to determine active put/get timesteps.
CALL mpi_bcast( pdata, nitems, MPI_WORKING_PRECISION, localRoot, mpi_comm_oce, ierr )
#endif
!
END SUBROUTINE cpl_rcv_1d
INTEGER FUNCTION cpl_freq( cdfieldname )
!!---------------------------------------------------------------------
......@@ -500,6 +812,7 @@ CONTAINS
!!----------------------------------------------------------------------
!
DEALLOCATE( exfld )
DEALLOCATE( exfld_ext )
IF(nstop == 0) THEN
CALL oasis_terminate( nerror )
ELSE
......@@ -543,7 +856,7 @@ CONTAINS
SUBROUTINE oasis_def_var(k1,cd1,k2,k3,k4,k5,k6,k7)
CHARACTER(*), INTENT(in ) :: cd1
INTEGER , INTENT(in ) :: k2,k3(2),k4,k5(2,2),k6
INTEGER , INTENT(in ) :: k2,k3(2),k4,k5(*),k6
INTEGER , INTENT( out) :: k1,k7
k1 = -1 ; k7 = -1
WRITE(numout,*) 'oasis_def_var: Error you sould not be there...', cd1
......@@ -563,13 +876,21 @@ CONTAINS
WRITE(numout,*) 'oasis_put: Error you sould not be there...'
END SUBROUTINE oasis_put
SUBROUTINE oasis_get(k1,k2,p1,k3)
SUBROUTINE oasis_get_2d(k1,k2,p1,k3)
REAL(wp), DIMENSION(:,:), INTENT( out) :: p1
INTEGER , INTENT(in ) :: k1,k2
INTEGER , INTENT( out) :: k3
p1(1,1) = -1. ; k3 = -1
WRITE(numout,*) 'oasis_get: Error you sould not be there...'
END SUBROUTINE oasis_get
WRITE(numout,*) 'oasis_get_2d: Error you sould not be there...'
END SUBROUTINE oasis_get_2d
SUBROUTINE oasis_get_1d(k1,k2,p1,k3)
REAL(wp), DIMENSION(:) , INTENT( out) :: p1
INTEGER , INTENT(in ) :: k1,k2
INTEGER , INTENT( out) :: k3
p1(1) = -1. ; k3 = -1
WRITE(numout,*) 'oasis_get_1d: Error you sould not be there...'
END SUBROUTINE oasis_get_1d
SUBROUTINE oasis_get_freqs(k1,k5,k2,k3,k4)
INTEGER , INTENT(in ) :: k1,k2
......
src/OCE/SBC/cpl_rnf_1d.F90 0 → 100644
View file @ 672ff863
MODULE cpl_rnf_1d
!!======================================================================
!! *** MODULE cpl_rnf_1d ***
!! Ocean forcing: River runoff passed from the atmosphere using
!! 1D array. One value per river.
!!=====================================================================
!! History : ?.? ! 2018-01 (D. Copsey) Initial setup
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! cpl_rnf_1d_init : runoffs initialisation
!!----------------------------------------------------------------------
#if defined key_oasis3
USE mod_oasis ! OASIS3-MCT module
#endif
USE timing ! Timing
USE in_out_manager ! I/O units
USE lib_mpp ! MPP library
USE iom
USE dom_oce ! Domain sizes (for grid box area e1e2t)
USE sbc_oce ! Surface boundary condition: ocean fields
USE lib_fortran, ONLY: DDPDD
IMPLICIT NONE
PRIVATE
PUBLIC cpl_rnf_1d_init ! routine called in nemo_init
PUBLIC cpl_rnf_1d_to_2d ! routine called in sbccpl.F90
TYPE, PUBLIC :: RIVERS_DATA !: Storage for river outflow data
INTEGER, ALLOCATABLE, DIMENSION(:,:) :: river_number !: River outflow number
REAL(wp), ALLOCATABLE, DIMENSION(:) :: river_area ! 1D array listing areas of each river outflow (m2)
COMPLEX(wp), ALLOCATABLE, DIMENSION(:) :: river_area_c ! Comlex version of river_area for use in bit reproducible sums (m2)
END TYPE RIVERS_DATA
! DB: see https://forge.ipsl.jussieu.fr/nemo/ticket/1865
#if defined key_agrif
TYPE(RIVERS_DATA), PUBLIC :: rivers !: River data
#else
TYPE(RIVERS_DATA), PUBLIC, TARGET :: rivers !: River data
#endif
INTEGER, PUBLIC :: nn_cpl_river ! Maximum number of rivers being passed through the coupler
INTEGER, PUBLIC :: runoff_id ! OASIS coupling id used in oasis_get command
LOGICAL :: ln_print_river_info ! Diagnostic prints of river coupling information
CONTAINS
SUBROUTINE cpl_rnf_1d_init
!!----------------------------------------------------------------------
!! *** SUBROUTINE cpl_rnf_1d_init ***
!!
!! ** Purpose : - Read in file for river outflow numbers.
!! Calculate 2D area of river outflow points.
!! Called from nemo_init (nemogcm.F90).
!!
!!----------------------------------------------------------------------
!! namelist variables
!!-------------------
CHARACTER(len=200) :: file_riv_number !: Filename for river numbers
INTEGER :: ios ! Local integer output status for namelist read
INTEGER :: inum
INTEGER :: ii, jj, rr !: Loop indices
INTEGER :: max_river
REAL(wp), DIMENSION(jpi,jpj) :: river_number ! 2D array containing the river outflow numbers
NAMELIST/nam_cpl_rnf_1d/file_riv_number, nn_cpl_river, ln_print_river_info
!!----------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('cpl_rnf_1d_init')
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'cpl_rnf_1d_init : initialization of river runoff coupling'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
IF(lwp) CALL flush(numout)
!! REWIND(numnam_cfg)
! Read the namelist
READ ( numnam_ref, nam_cpl_rnf_1d, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_cpl_rnf_1d in reference namelist' )
READ ( numnam_cfg, nam_cpl_rnf_1d, IOSTAT = ios, ERR = 902 )
902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_cpl_rnf_1d in configuration namelist' )
IF(lwm) WRITE ( numond, nam_cpl_rnf_1d )
! ! Parameter control and print
IF(lwp) WRITE(numout,*) ' '
IF(lwp) WRITE(numout,*) ' Namelist nam_cpl_rnf_1d : Coupled runoff using 1D array'
IF(lwp) WRITE(numout,*) ' Input file that contains river numbers = ',file_riv_number
IF(lwp) WRITE(numout,*) ' Maximum number of rivers to couple = ',nn_cpl_river
IF(lwp) WRITE(numout,*) ' Print river information = ',ln_print_river_info
IF(lwp) WRITE(numout,*) ' '
IF(lwp) CALL flush(numout)
! Assign space for river numbers
ALLOCATE( rivers%river_number( jpi, jpj ) )
! Read the river numbers from netcdf file
CALL iom_open (file_riv_number , inum )
CALL iom_get ( inum, jpdom_global, 'river_number', river_number, cd_type = 'T', psgn = 1._wp)
CALL iom_close( inum )
! Convert from a real array to an integer array
max_river=0
DO ii = 1, jpi
DO jj = 1, jpj
rivers%river_number(ii,jj) = INT(river_number(ii,jj))
IF ( rivers%river_number(ii,jj) > max_river ) THEN
max_river = rivers%river_number(ii,jj)
END IF
END DO
END DO
! Print out the largest river number
IF ( ln_print_river_info .AND. lwp) THEN
WRITE(numout,*) 'Maximum river number in input file = ',max_river
CALL flush(numout)
END IF
! Get the area of each river outflow
ALLOCATE( rivers%river_area( nn_cpl_river ) )
ALLOCATE( rivers%river_area_c( nn_cpl_river ) )
rivers%river_area_c(:) = CMPLX( 0.e0, 0.e0, wp )
DO ii = Nis0, Nie0
DO jj = Njs0, Nje0
IF ( tmask_i(ii,jj) > 0.5 ) THEN ! This makes sure we are not at a duplicated point (at north fold or east-west cyclic point)
IF ( rivers%river_number(ii,jj) > 0 .AND. rivers%river_number(ii,jj) <= nn_cpl_river ) THEN
! Add the area of each grid box (e1e2t) into river_area_c using DDPDD which should maintain bit reproducibility (needs to be checked)
CALL DDPDD( CMPLX( e1e2t(ii,jj), 0.e0, wp ), rivers%river_area_c(rivers%river_number(ii,jj) ) )
END IF
END IF
END DO
END DO
! Use mpp_sum to add together river areas on other processors
CALL mpp_sum( 'cpl_rnf_1d', rivers%river_area_c )
! Convert from complex number to real
rivers%river_area(:) = REAL(rivers%river_area_c(:),wp)
IF ( ln_print_river_info .AND. lwp) THEN
WRITE(numout,*) 'Area of rivers 1 to 10 are ',rivers%river_area(1:10)
WRITE(numout,*) 'Maximum river area = ',MAXVAL(rivers%river_area)
WRITE(numout,*) 'Minimum river area = ',MINVAL(rivers%river_area)
WRITE(numout,*) 'Sum of river areas = ',SUM(rivers%river_area)
CALL flush(numout)
END IF
IF ( MINVAL(rivers%river_area) <= 0 ) THEN
WRITE(numout,*) 'ERROR: There is at least one river with a river outflow area of zero. Please check your file_riv_number file'
WRITE(numout,*) 'that all the allocated river numbers are at ocean points as defined by the bathymetry file with no river'
WRITE(numout,*) 'numbers within the north fold or wraparound points.'
DO rr = 1,nn_cpl_river
IF ( rivers%river_area(rr) <= 0 ) THEN
WRITE(numout,*) 'The first river with this problem is river number ',rr
CALL ctl_stop ( 'STOP', 'ERROR: There is at least one river with a river outflow area of zero. See stdout.')
END IF
END DO
END IF
END SUBROUTINE cpl_rnf_1d_init
SUBROUTINE cpl_rnf_1d_to_2d( runoff_1d )
!!----------------------------------------------------------------------
!! *** SUBROUTINE cpl_rnf_1d_to_2d ***
!!
!! ** Purpose : - Convert river outflow from 1D array (passed from the
!! atmosphere) to the 2D NEMO runoff field.
!! Called from sbc_cpl_ice_flx (sbccpl.F90).
!!
!!----------------------------------------------------------------------
REAL , INTENT(in ) :: runoff_1d(nn_cpl_river) ! River runoff. One value per river.
INTEGER :: ii, jj, rr ! Loop indices
LOGICAL :: found_nan
! Convert the 1D total runoff per river to 2D runoff flux by
! dividing by the area of each runoff zone.
DO ii = 1, jpi
DO jj = 1, jpj
IF ( rivers%river_number(ii,jj) > 0 .AND. rivers%river_number(ii,jj) <= nn_cpl_river .AND. tmask_i(ii,jj) > 0.5 ) THEN
rnf(ii,jj) = runoff_1d(rivers%river_number(ii,jj)) / rivers%river_area(rivers%river_number(ii,jj))
ELSE
rnf(ii,jj) = 0.0
END IF
END DO
END DO
IF ( ln_print_river_info ) THEN
WRITE(numout,*) 'SUM of river runoff on 1D array = ',SUM(runoff_1d)
WRITE(numout,*) 'SUM of river runoff on 2D array = ',SUM(rnf)
found_nan = .false.
DO rr = 1, nn_cpl_river
IF (rr <= 10) THEN
WRITE(numout,*) 'River number ',rr,' has total runoff=',runoff_1d(rr),' and area=',rivers%river_area(rr),' making runoff flux=',runoff_1d(rr)/rivers%river_area(rr)
END IF
IF (runoff_1d(rr) /= runoff_1d(rr)) THEN
IF (.NOT. found_nan) THEN
WRITE(numout,*) 'WARNING: NAN found at river number ',rr
found_nan = .true.
END IF
END IF
END DO
END IF
END SUBROUTINE cpl_rnf_1d_to_2d
END MODULE cpl_rnf_1d
src/OCE/SBC/geo2ocean.F90
View file @ 672ff863
......@@ -25,6 +25,7 @@ MODULE geo2ocean
IMPLICIT NONE
PRIVATE
PUBLIC repcmo ! called in sbccpl
PUBLIC rot_rep ! called in sbccpl, fldread, and cyclone
PUBLIC geo2oce ! called in sbccpl
PUBLIC oce2geo ! called in sbccpl
......@@ -50,6 +51,48 @@ MODULE geo2ocean
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE repcmo ( pxu1, pyu1, pxv1, pyv1, &
px2 , py2 , kchoix )
!!----------------------------------------------------------------------
!! *** ROUTINE repcmo ***
!!
!! ** Purpose : Change vector componantes from a geographic grid to a
!! stretched coordinates grid.
!!
!! ** Method : Initialization of arrays at the first call.
!!
!! ** Action : - px2 : first componante (defined at u point)
!! - py2 : second componante (defined at v point)
!!----------------------------------------------------------------------
REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pxu1, pyu1 ! geographic vector componantes at u-point
REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pxv1, pyv1 ! geographic vector componantes at v-point
REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: px2 ! i-componante (defined at u-point)
REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: py2 ! j-componante (defined at v-point)
!!----------------------------------------------------------------------
INTEGER, INTENT( IN ) :: &
kchoix ! type of transformation
! = 1 change from geographic to model grid.
! =-1 change from model to geographic grid
!!----------------------------------------------------------------------
SELECT CASE (kchoix)
CASE ( 1)
! Change from geographic to stretched coordinate
! ----------------------------------------------
CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 )
CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 )
CASE (-1)
! Change from stretched to geographic coordinate
! ----------------------------------------------
CALL rot_rep( pxu1, pyu1, 'U', 'ij->e',px2 )
CALL rot_rep( pxv1, pyv1, 'V', 'ij->n',py2 )
END SELECT
END SUBROUTINE repcmo
SUBROUTINE rot_rep ( pxin, pyin, cd_type, cdtodo, prot )
!!----------------------------------------------------------------------
!! *** ROUTINE rot_rep ***
......
src/OCE/SBC/sbc_ice.F90
View file @ 672ff863
......@@ -94,6 +94,7 @@ MODULE sbc_ice
! already defined in ice.F90 for SI3
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: h_i, h_s
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_last_couple !: Sea ice fraction on categories at the last coupling point
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tatm_ice !: air temperature [K]
#endif
......
src/OCE/SBC/sbc_oce.F90
View file @ 672ff863
......@@ -131,6 +131,8 @@ MODULE sbc_oce
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: atm_co2 !: atmospheric pCO2 [ppm]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: cloud_fra !: cloud cover (fraction of cloud in a gridcell) [-]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: greenland_icesheet_mask
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: antarctica_icesheet_mask
!!---------------------------------------------------------------------
!! ABL Vertical Domain size
......@@ -153,6 +155,33 @@ MODULE sbc_oce
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3t_m !: mean (nn_fsbc time-step) sea surface layer thickness [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_m !: mean (nn_fsbc time-step) fraction of solar net radiation absorbed in the 1st T level [-]
!!----------------------------------------------------------------------
!! Surface scalars of total ice sheet mass for Greenland and Antarctica,
!! passed from atmosphere to be converted to dvol and hence a freshwater
!! flux by using old values. New values are saved in the dump, to become
!! old values next coupling timestep. Freshwater fluxes split between
!! sub iceshelf melting and iceberg calving, scalled to flux per second
!!----------------------------------------------------------------------
REAL(wp), PUBLIC :: greenland_icesheet_mass, greenland_icesheet_mass_rate_of_change, greenland_icesheet_timelapsed
REAL(wp), PUBLIC :: antarctica_icesheet_mass, antarctica_icesheet_mass_rate_of_change, antarctica_icesheet_timelapsed
! sbccpl namelist parameters associated with icesheet freshwater input code. Included here rather than in sbccpl.F90 to
! avoid circular dependencies.
INTEGER, PUBLIC :: nn_coupled_iceshelf_fluxes ! =0 : total freshwater input from iceberg calving and ice shelf basal melting
! taken from climatologies used (no action in coupling routines).
! =1 : use rate of change of mass of Greenland and Antarctic icesheets to set the
! combined magnitude of the iceberg calving and iceshelf melting freshwater fluxes.
! =2 : specify constant freshwater inputs in this namelist to set the combined
! magnitude of iceberg calving and iceshelf melting freshwater fluxes.
LOGICAL, PUBLIC :: ln_iceshelf_init_atmos ! If true force ocean to initialise iceshelf masses from atmospheric values rather
! than values in ocean restart (applicable if nn_coupled_iceshelf_fluxes=1).
REAL(wp), PUBLIC :: rn_greenland_total_fw_flux ! Constant total rate of freshwater input (kg/s) for Greenland (if nn_coupled_iceshelf_fluxes=2)
REAL(wp), PUBLIC :: rn_greenland_calving_fraction ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
REAL(wp), PUBLIC :: rn_antarctica_total_fw_flux ! Constant total rate of freshwater input (kg/s) for Antarctica (if nn_coupled_iceshelf_fluxes=2)
REAL(wp), PUBLIC :: rn_antarctica_calving_fraction ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
REAL(wp), PUBLIC :: rn_iceshelf_fluxes_tolerance ! Absolute tolerance for detecting differences in icesheet masses.
!!----------------------------------------------------------------------
!! Surface atmospheric fields
!!----------------------------------------------------------------------
......@@ -193,6 +222,9 @@ CONTAINS
& atm_co2(jpi,jpj) , tsk_m(jpi,jpj) , cloud_fra(jpi,jpj), &
& ssu_m (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) , &
& ssv_m (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) )
ALLOCATE( greenland_icesheet_mask(jpi,jpj) , antarctica_icesheet_mask(jpi,jpj) )
!
ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) )
!
......
src/OCE/SBC/sbcblk.F90
View file @ 672ff863
......@@ -119,6 +119,7 @@ MODULE sbcblk
!
REAL(dp) :: rn_pfac ! multiplication factor for precipitation
REAL(wp), PUBLIC :: rn_efac ! multiplication factor for evaporation
REAL(wp) :: rn_vfac ! multiplication factor for ice/ocean velocity in the calculation of wind stress !
REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements
REAL(wp) :: rn_zu ! z(u) : height of wind measurements
!
......@@ -220,7 +221,7 @@ CONTAINS
INTEGER :: ipka ! number of levels in the atmospheric variable
NAMELIST/namsbc_blk/ ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, ln_ANDREAS, & ! bulk algorithm
& rn_zqt, rn_zu, nn_iter_algo, ln_skin_cs, ln_skin_wl, &
& rn_pfac, rn_efac, &
& rn_pfac, rn_efac, rn_vfac, &
& ln_crt_fbk, rn_stau_a, rn_stau_b, & ! current feedback
& ln_humi_sph, ln_humi_dpt, ln_humi_rlh, ln_tair_pot, &
& ln_Cx_ice_cst, rn_Cd_i, rn_Ce_i, rn_Ch_i, &
......@@ -414,6 +415,7 @@ CONTAINS
WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu
WRITE(numout,*) ' factor applied on precipitation (total & snow) rn_pfac = ', rn_pfac
WRITE(numout,*) ' factor applied on evaporation rn_efac = ', rn_efac
WRITE(numout,*) ' factor applied on ocean/ice velocity rn_vfac = ', rn_vfac
WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))'
WRITE(numout,*) ' use surface current feedback on wind stress ln_crt_fbk = ', ln_crt_fbk
IF(ln_crt_fbk) THEN
......@@ -655,9 +657,7 @@ CONTAINS
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zztmp ! local variable
REAL(wp) :: zstmax, zstau
#if defined key_cyclone
REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point
#endif
REAL(wp), DIMENSION(jpi,jpj) :: ztau_i, ztau_j ! wind stress components at T-point
REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s]
REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean
......@@ -683,9 +683,10 @@ CONTAINS
! ----------------------------------------------------------------------------- !
! ... components ( U10m - U_oce ) at T-point (unmasked)
#if defined key_cyclone
zwnd_i(:,:) = 0._wp
zwnd_j(:,:) = 0._wp
#if defined key_cyclone
CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj)
......@@ -695,9 +696,15 @@ CONTAINS
END_2D
#else
! ... scalar wind module at T-point (not masked)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
wndm(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
DO_2D( 0, 0, 0, 0 )
zwnd_i(ji,jj) = ( pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) )
zwnd_j(ji,jj) = ( pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) )
END_2D
CALL lbc_lnk( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. )
! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) &
& + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1)
#endif
! ----------------------------------------------------------------------------- !
! I Solar FLUX !
......@@ -827,8 +834,13 @@ CONTAINS
ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
#else
ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
IF ( rn_vfac > 0._wp ) THEN
ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
ELSE
ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
ENDIF
#endif
ELSE
ztau_i(ji,jj) = 0._wp
......@@ -1017,7 +1029,8 @@ CONTAINS
!
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zootm_su ! sea-ice surface mean temperature
REAL(wp) :: zztmp1, zztmp2 ! temporary scalars
REAL(wp) :: zwndi_t , zwndj_t
REAL(wp) :: zztmp0,zztmp1, zztmp2 ! temporary scalars
REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zsipt ! temporary array
!!---------------------------------------------------------------------
!
......@@ -1025,9 +1038,12 @@ CONTAINS
! Wind module relative to the moving ice ( U10m - U_ice ) !
! ------------------------------------------------------------ !
! C-grid ice dynamics : U & V-points (same as ocean)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
DO_2D( 0, 0, 0, 0 )
zwndi_t = ( pwndi(ji,jj) - rn_vfac * 0.5 * ( puice(ji-1,jj ) + puice(ji,jj) ) )
zwndj_t = ( pwndj(ji,jj) - rn_vfac * 0.5 * ( pvice(ji ,jj-1) + pvice(ji,jj) ) )
wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t )* tmask(ji,jj,1)
END_2D
CALL lbc_lnk( 'sbcblk', wndm_ice, 'T', 1. )
!
! potential sea-ice surface temperature [K]
zsipt(:,:) = theta_exner( ptsui(:,:), pslp(:,:) )
......@@ -1074,10 +1090,11 @@ CONTAINS
! Wind stress relative to nonmoving ice ( U10m ) !
! ---------------------------------------------------- !
! supress moving ice in wind stress computation as we don't know how to do it properly...
zztmp0 = rn_vfac * 0.5_wp
DO_2D( 0, 1, 0, 1 ) ! at T point
zztmp1 = rhoa(ji,jj) * Cd_ice(ji,jj) * wndm_ice(ji,jj)
putaui(ji,jj) = zztmp1 * pwndi(ji,jj)
pvtaui(ji,jj) = zztmp1 * pwndj(ji,jj)
putaui(ji,jj) = zztmp1 * ( pwndi(ji,jj) - zztmp0 * ( puice(ji-1,jj ) + puice(ji,jj) ) )
pvtaui(ji,jj) = zztmp1 * ( pwndj(ji,jj) - zztmp0 * ( pvice(ji ,jj-1) + pvice(ji,jj) ) )
END_2D
!#LB: saving the module, and x-y components, of the ai wind-stress at T-points: NOT weighted by the ice concentration !!!
......
src/OCE/SBC/sbccpl.F90
View file @ 672ff863
......@@ -33,10 +33,18 @@ MODULE sbccpl
#endif
USE cpl_oasis3 ! OASIS3 coupling
USE geo2ocean !
USE oce , ONLY : ts, uu, vv, ssh, fraqsr_1lev
USE oce
#if defined key_medusa
USE oce , ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl, &
PCO2a_in_cpl, Dust_in_cpl
#endif
USE ocealb !
USE eosbn2 !
USE sbcrnf , ONLY : l_rnfcpl
USE cpl_rnf_1d, ONLY: nn_cpl_river, cpl_rnf_1d_init, cpl_rnf_1d_to_2d ! Variables used in 1D river outflow
#if defined key_medusa
USE par_trc , ONLY : ln_medusa
#endif
#if defined key_cice
USE ice_domain_size, only: ncat
#endif
......@@ -48,6 +56,7 @@ MODULE sbccpl
USE iom ! NetCDF library
USE lib_mpp ! distribued memory computing library
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE timing
#if defined key_oasis3
USE mod_oasis, ONLY : OASIS_Sent, OASIS_ToRest, OASIS_SentOut, OASIS_ToRestOut
......@@ -129,8 +138,17 @@ MODULE sbccpl
INTEGER, PARAMETER :: jpr_icb = 61
INTEGER, PARAMETER :: jpr_ts_ice = 62 ! Sea ice surface temp
!!INTEGER, PARAMETER :: jpr_qtrice = 63 ! Transmitted solar thru sea-ice
INTEGER, PARAMETER :: jprcv = 62 ! total number of fields received
INTEGER, PARAMETER :: jpr_grnm = 63 ! Greenland ice mass
INTEGER, PARAMETER :: jpr_antm = 64 ! Antarctic ice mass
INTEGER, PARAMETER :: jpr_rnf_1d = 65 ! 1D river runoff
INTEGER, PARAMETER :: jpr_qtr = 66 ! Transmitted solar
#if defined key_medusa
INTEGER, PARAMETER :: jpr_atm_pco2 = 67 ! Incoming atm pCO2 flux
INTEGER, PARAMETER :: jpr_atm_dust = 68 ! Incoming atm aggregate dust
INTEGER, PARAMETER :: jprcv = 69 ! total number of fields received
#else
INTEGER, PARAMETER :: jprcv = 66 ! total number of fields received
#endif
INTEGER, PARAMETER :: jps_fice = 1 ! ice fraction sent to the atmosphere
INTEGER, PARAMETER :: jps_toce = 2 ! ocean temperature
......@@ -170,17 +188,24 @@ MODULE sbccpl
INTEGER, PARAMETER :: jps_kice = 36 ! sea ice effective conductivity
INTEGER, PARAMETER :: jps_sstfrz = 37 ! sea surface freezing temperature
INTEGER, PARAMETER :: jps_ttilyr = 38 ! sea ice top layer temp
INTEGER, PARAMETER :: jpsnd = 38 ! total number of fields sent
#if ! defined key_oasis3
! Dummy variables to enable compilation when oasis3 is not being used
INTEGER :: OASIS_Sent = -1
INTEGER :: OASIS_SentOut = -1
INTEGER :: OASIS_ToRest = -1
INTEGER :: OASIS_ToRestOut = -1
#if defined key_medusa
INTEGER, PARAMETER :: jps_bio_co2 = 39 ! MEDUSA air-sea CO2 flux
INTEGER, PARAMETER :: jps_bio_dms = 40 ! MEDUSA DMS surface concentration
INTEGER, PARAMETER :: jps_bio_chloro = 41 ! MEDUSA chlorophyll surface concentration
INTEGER, PARAMETER :: jpsnd = 41 ! total number of fields sent
#else
INTEGER, PARAMETER :: jpsnd = 38 ! total number of fields sent
#endif
#if ! defined key_oasis3
! Dummy variables to enable compilation when oasis3 is not being used
INTEGER :: OASIS_Sent = -1
INTEGER :: OASIS_SentOut = -1
INTEGER :: OASIS_ToRest = -1
INTEGER :: OASIS_ToRestOut = -1
#endif
! !!** namelist namsbc_cpl **
TYPE :: FLD_C !
CHARACTER(len = 32) :: cldes ! desciption of the coupling strategy
......@@ -193,20 +218,34 @@ MODULE sbccpl
TYPE(FLD_C) :: sn_snd_temp , sn_snd_alb , sn_snd_thick, sn_snd_crt , sn_snd_co2, &
& sn_snd_thick1, sn_snd_cond, sn_snd_mpnd , sn_snd_sstfrz, sn_snd_ttilyr
! ! Received from the atmosphere
TYPE(FLD_C) :: sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, &
#if defined key_medusa
TYPE(FLD_C) :: sn_snd_bio_co2, sn_snd_bio_dms, sn_snd_bio_chloro
#endif
TYPE(FLD_C) :: sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_tauw, sn_rcv_dqnsdt, sn_rcv_qsr, &
& sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf, sn_rcv_ts_ice
TYPE(FLD_C) :: sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf
! ! Send to waves
TYPE(FLD_C) :: sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev
! ! Received from waves
TYPE(FLD_C) :: sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf, &
sn_rcv_grnm, sn_rcv_antm
#if defined key_medusa
TYPE(FLD_C) :: sn_rcv_atm_pco2, sn_rcv_atm_dust
#endif
! Send to waves
TYPE(FLD_C) :: sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev
TYPE(FLD_C) :: sn_rcv_hsig, sn_rcv_phioc, sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper, sn_rcv_wnum, &
& sn_rcv_wstrf, sn_rcv_wdrag, sn_rcv_charn, sn_rcv_taw, sn_rcv_bhd, sn_rcv_tusd, sn_rcv_tvsd
! Received from waves
TYPE(FLD_C) :: sn_rcv_tauwoc, sn_rcv_wfreq
! Transmitted solar
TYPE(FLD_C) :: sn_rcv_qtr
! ! Other namelist parameters
!! TYPE(FLD_C) :: sn_rcv_qtrice
INTEGER :: nn_cplmodel ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
LOGICAL :: ln_usecplmask ! use a coupling mask file to merge data received from several models
! -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
LOGICAL :: ln_scale_ice_flux ! use ice fluxes that are already "ice weighted" ( i.e. multiplied ice concentration)
LOGICAL :: ln_couple_ocean_evap ! Do we couple total (ocean+sea ice) evaporation (FALSE)
! or ocean only evaporation (TRUE)
TYPE :: DYNARR
REAL(wp), POINTER, DIMENSION(:,:,:) :: z3
......@@ -215,9 +254,7 @@ MODULE sbccpl
TYPE( DYNARR ), SAVE, DIMENSION(jprcv) :: frcv ! all fields recieved from the atmosphere
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: alb_oce_mix ! ocean albedo sent to atmosphere (mix clear/overcast sky)
#if defined key_si3 || defined key_cice
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_last_couple !: Ice fractional area at last coupling time
#endif
INTEGER , ALLOCATABLE, SAVE, DIMENSION(:) :: nrcvinfo ! OASIS info argument
......@@ -275,22 +312,38 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj) :: zacs, zaos
!!
NAMELIST/namsbc_cpl/ nn_cplmodel , ln_usecplmask, nn_cats_cpl , ln_scale_ice_flux, &
& ln_couple_ocean_evap, &
& sn_snd_temp , sn_snd_alb , sn_snd_thick, sn_snd_crt , sn_snd_co2 , &
& sn_snd_ttilyr, sn_snd_cond , sn_snd_mpnd , sn_snd_sstfrz, sn_snd_thick1, &
& sn_snd_ifrac , sn_snd_crtw , sn_snd_wlev , sn_rcv_hsig , sn_rcv_phioc , &
& sn_rcv_w10m , sn_rcv_taumod, sn_rcv_tau , sn_rcv_dqnsdt, sn_rcv_qsr , &
& sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum , sn_rcv_wstrf , &
& sn_rcv_charn , sn_rcv_taw , sn_rcv_bhd , sn_rcv_tusd , sn_rcv_tvsd, &
sn_rcv_qtr , &
& sn_rcv_wdrag , sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf , sn_rcv_cal , &
& sn_rcv_iceflx, sn_rcv_co2 , sn_rcv_icb , sn_rcv_isf , sn_rcv_ts_ice, & !!, sn_rcv_qtrice
& sn_rcv_mslp
& sn_rcv_mslp, &
sn_rcv_grnm , sn_rcv_antm , &
& nn_coupled_iceshelf_fluxes , ln_iceshelf_init_atmos , &
& rn_greenland_total_fw_flux , rn_greenland_calving_fraction , &
& rn_antarctica_total_fw_flux , rn_antarctica_calving_fraction , &
#if defined key_medusa
& rn_iceshelf_fluxes_tolerance, &
! Add MEDUSA related fields to namelist
sn_snd_bio_co2 , sn_snd_bio_dms, sn_snd_bio_chloro, &
& sn_rcv_atm_pco2, sn_rcv_atm_dust
#else
& rn_iceshelf_fluxes_tolerance
#endif
!!---------------------------------------------------------------------
!
! ================================ !
! Namelist informations !
! ================================ !
!
!
IF (ln_timing) CALL timing_start('sbc_cpl_init')
READ ( numnam_ref, namsbc_cpl, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in reference namelist' )
!
......@@ -307,6 +360,7 @@ CONTAINS
WRITE(numout,*)' nn_cplmodel = ', nn_cplmodel
WRITE(numout,*)' ln_usecplmask = ', ln_usecplmask
WRITE(numout,*)' ln_scale_ice_flux = ', ln_scale_ice_flux
WRITE(numout,*)' ln_couple_ocean_evap = ', ln_couple_ocean_evap
WRITE(numout,*)' nn_cats_cpl = ', nn_cats_cpl
WRITE(numout,*)' received fields (mutiple ice categogies)'
WRITE(numout,*)' 10m wind module = ', TRIM(sn_rcv_w10m%cldes ), ' (', TRIM(sn_rcv_w10m%clcat ), ')'
......@@ -321,10 +375,12 @@ CONTAINS
WRITE(numout,*)' freshwater budget = ', TRIM(sn_rcv_emp%cldes ), ' (', TRIM(sn_rcv_emp%clcat ), ')'
WRITE(numout,*)' runoffs = ', TRIM(sn_rcv_rnf%cldes ), ' (', TRIM(sn_rcv_rnf%clcat ), ')'
WRITE(numout,*)' calving = ', TRIM(sn_rcv_cal%cldes ), ' (', TRIM(sn_rcv_cal%clcat ), ')'
WRITE(numout,*)' Greenland ice mass = ', TRIM(sn_rcv_grnm%cldes ), ' (', TRIM(sn_rcv_grnm%clcat ), ')'
WRITE(numout,*)' Antarctica ice mass = ', TRIM(sn_rcv_antm%cldes ), ' (', TRIM(sn_rcv_antm%clcat ), ')'
WRITE(numout,*)' iceberg = ', TRIM(sn_rcv_icb%cldes ), ' (', TRIM(sn_rcv_icb%clcat ), ')'
WRITE(numout,*)' ice shelf = ', TRIM(sn_rcv_isf%cldes ), ' (', TRIM(sn_rcv_isf%clcat ), ')'
WRITE(numout,*)' sea ice heat fluxes = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
!! WRITE(numout,*)' transmitted solar thru sea-ice = ', TRIM(sn_rcv_qtrice%cldes), ' (', TRIM(sn_rcv_qtrice%clcat), ')'
WRITE(numout,*)' transmitted solar = ', TRIM(sn_rcv_qtr%cldes ), ' (', TRIM(sn_rcv_qtr%clcat ), ')'
WRITE(numout,*)' atm co2 = ', TRIM(sn_rcv_co2%cldes ), ' (', TRIM(sn_rcv_co2%clcat ), ')'
WRITE(numout,*)' Sea ice surface skin temperature= ', TRIM(sn_rcv_ts_ice%cldes), ' (', TRIM(sn_rcv_ts_ice%clcat), ')'
WRITE(numout,*)' surface waves:'
......@@ -347,6 +403,11 @@ CONTAINS
WRITE(numout,*)' - referential = ', sn_snd_crt%clvref
WRITE(numout,*)' - orientation = ', sn_snd_crt%clvor
WRITE(numout,*)' - mesh = ', sn_snd_crt%clvgrd
#if defined key_medusa
WRITE(numout,*)' bio co2 flux = ', TRIM(sn_snd_bio_co2%cldes),' (', TRIM(sn_snd_bio_co2%clcat), ')'
WRITE(numout,*)' bio dms flux = ', TRIM(sn_snd_bio_dms%cldes),' (', TRIM(sn_snd_bio_dms%clcat), ')'
WRITE(numout,*)' bio dms chlorophyll = ', TRIM(sn_snd_bio_chloro%cldes), ' (', TRIM(sn_snd_bio_chloro%clcat), ')'
#endif
WRITE(numout,*)' oce co2 flux = ', TRIM(sn_snd_co2%cldes ), ' (', TRIM(sn_snd_co2%clcat ), ')'
WRITE(numout,*)' ice effective conductivity = ', TRIM(sn_snd_cond%cldes ), ' (', TRIM(sn_snd_cond%clcat ), ')'
WRITE(numout,*)' meltponds fraction and depth = ', TRIM(sn_snd_mpnd%cldes ), ' (', TRIM(sn_snd_mpnd%clcat ), ')'
......@@ -357,6 +418,17 @@ CONTAINS
WRITE(numout,*)' - referential = ', sn_snd_crtw%clvref
WRITE(numout,*)' - orientation = ', sn_snd_crtw%clvor
WRITE(numout,*)' - mesh = ', sn_snd_crtw%clvgrd
#if defined key_medusa
WRITE(numout,*)' atm pco2 = ', TRIM(sn_rcv_atm_pco2%cldes),'(', TRIM(sn_rcv_atm_pco2%clcat ), ')'
WRITE(numout,*)' atm dust = ', TRIM(sn_rcv_atm_dust%cldes),'(', TRIM(sn_rcv_atm_dust%clcat),')'
#endif
WRITE(numout,*)' nn_coupled_iceshelf_fluxes = ', nn_coupled_iceshelf_fluxes
WRITE(numout,*)' ln_iceshelf_init_atmos = ', ln_iceshelf_init_atmos
WRITE(numout,*)' rn_greenland_total_fw_flux = ', rn_greenland_total_fw_flux
WRITE(numout,*)' rn_antarctica_total_fw_flux = ', rn_antarctica_total_fw_flux
WRITE(numout,*)' rn_greenland_calving_fraction = ', rn_greenland_calving_fraction
WRITE(numout,*)' rn_antarctica_calving_fraction = ', rn_antarctica_calving_fraction
WRITE(numout,*)' rn_iceshelf_fluxes_tolerance = ', rn_iceshelf_fluxes_tolerance
ENDIF
IF( lwp .AND. ln_wave) THEN ! control print
WRITE(numout,*)' surface waves:'
......@@ -391,7 +463,11 @@ CONTAINS
! define the north fold type of lbc (srcv(:)%nsgn)
! default definitions of srcv
srcv(:)%laction = .FALSE. ; srcv(:)%clgrid = 'T' ; srcv(:)%nsgn = 1. ; srcv(:)%nct = 1
srcv(:)%laction = .FALSE.
srcv(:)%clgrid = 'T'
srcv(:)%nsgn = 1.
srcv(:)%nct = 1
srcv(:)%dimensions = 2
! ! ------------------------- !
! ! ice and ocean wind stress !
......@@ -445,7 +521,11 @@ CONTAINS
srcv(jpr_otx1:jpr_otz1)%clgrid = 'U' ! oce components given at U-point
srcv(jpr_otx2:jpr_otz2)%clgrid = 'V' ! and V-point
srcv(jpr_itx1:jpr_itz1)%clgrid = 'F' ! ice components given at F-point
srcv(jpr_otx1:jpr_otz2)%laction = .TRUE. ! receive oce components on grid 1 & 2
!srcv(jpr_otx1:jpr_otz2)%laction = .TRUE. ! receive oce components on grid 1 & 2
! Currently needed for HadGEM3 - but shouldn't affect anyone else for the moment
srcv(jpr_otx1)%laction = .TRUE.
srcv(jpr_oty1)%laction = .TRUE.
!
srcv(jpr_itx1:jpr_itz1)%laction = .TRUE. ! receive ice components on grid 1 only
CASE( 'T,I' )
srcv(jpr_otx1:jpr_itz2)%clgrid = 'T' ! oce and ice components given at T-point
......@@ -510,15 +590,34 @@ CONTAINS
! ! Runoffs & Calving !
! ! ------------------------- !
srcv(jpr_rnf )%clname = 'O_Runoff'
IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN
srcv(jpr_rnf)%laction = .TRUE.
srcv(jpr_rnf_1d )%clname = 'ORunff1D'
IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' .OR. TRIM( sn_rcv_rnf%cldes ) == 'coupled1d' ) THEN
IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) srcv(jpr_rnf)%laction = .TRUE.
IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled1d' ) THEN
srcv(jpr_rnf_1d)%laction = .TRUE.
srcv(jpr_rnf_1d)%dimensions = 1 ! 1D field passed through coupler
END IF
l_rnfcpl = .TRUE. ! -> no need to read runoffs in sbcrnf
ln_rnf = nn_components /= jp_iam_sas ! -> force to go through sbcrnf if not sas
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' runoffs received from oasis -> force ln_rnf = ', ln_rnf
ENDIF
!
srcv(jpr_cal)%clname = 'OCalving' ; IF( TRIM( sn_rcv_cal%cldes) == 'coupled' ) srcv(jpr_cal)%laction = .TRUE.
srcv(jpr_cal )%clname = 'OCalving'
IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' ) srcv(jpr_cal)%laction = .TRUE.
srcv(jpr_grnm )%clname = 'OGrnmass'
IF( TRIM( sn_rcv_grnm%cldes ) == 'coupled' ) THEN
srcv(jpr_grnm)%laction = .TRUE.
srcv(jpr_grnm)%dimensions = 0 ! Scalar field
ENDIF
srcv(jpr_antm )%clname = 'OAntmass'
IF( TRIM( sn_rcv_antm%cldes ) == 'coupled' ) THEN
srcv(jpr_antm)%laction = .TRUE.
srcv(jpr_antm)%dimensions = 0 ! Scalar field
ENDIF
srcv(jpr_isf)%clname = 'OIcshelf' ; IF( TRIM( sn_rcv_isf%cldes) == 'coupled' ) srcv(jpr_isf)%laction = .TRUE.
srcv(jpr_icb)%clname = 'OIceberg' ; IF( TRIM( sn_rcv_icb%cldes) == 'coupled' ) srcv(jpr_icb)%laction = .TRUE.
......@@ -594,6 +693,22 @@ CONTAINS
IF(lwp) WRITE(numout,*) ' Atmospheric pco2 received from oasis '
IF(lwp) WRITE(numout,*)
ENDIF
#if defined key_medusa
! ! --------------------------------------- !
! ! Incoming CO2 and DUST fluxes for MEDUSA !
! ! --------------------------------------- !
srcv(jpr_atm_pco2)%clname = 'OATMPCO2'
IF (TRIM(sn_rcv_atm_pco2%cldes) == 'medusa') THEN
srcv(jpr_atm_pco2)%laction = .TRUE.
END IF
srcv(jpr_atm_dust)%clname = 'OATMDUST'
IF (TRIM(sn_rcv_atm_dust%cldes) == 'medusa') THEN
srcv(jpr_atm_dust)%laction = .TRUE.
END IF
#endif
!
! ! ------------------------- !
! ! Mean Sea Level Pressure !
......@@ -625,6 +740,21 @@ CONTAINS
!! ENDIF
!! srcv(jpr_qtrice)%laction = .TRUE.
!! ENDIF
! ! ------------------------- !
! ! transmitted solar !
! ! ------------------------- !
srcv(jpr_qtr )%clname = 'OQtr'
IF( TRIM(sn_rcv_qtr%cldes) == 'coupled' ) THEN
IF ( TRIM( sn_rcv_qtr%clcat ) == 'yes' ) THEN
srcv(jpr_qtr)%nct = nn_cats_cpl
ELSE
CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qtr%clcat should always be set to yes currently' )
ENDIF
srcv(jpr_qtr)%laction = .TRUE.
ENDIF
! ! ------------------------- !
! ! ice skin temperature !
! ! ------------------------- !
......@@ -799,24 +929,6 @@ CONTAINS
ENDIF
ENDIF
! =================================================== !
! Allocate all parts of frcv used for received fields !
! =================================================== !
DO jn = 1, jprcv
IF( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
END DO
! Allocate taum part of frcv which is used even when not received as coupling field
IF( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
! Allocate w10m part of frcv which is used even when not received as coupling field
IF( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) )
! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field
IF( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) )
IF( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) )
! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
IF( k_ice /= 0 ) THEN
IF( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) )
IF( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) )
ENDIF
! ================================ !
! Define the send interface !
......@@ -826,8 +938,12 @@ CONTAINS
! define the north fold type of lbc (ssnd(:)%nsgn)
! default definitions of nsnd
ssnd(:)%laction = .FALSE. ; ssnd(:)%clgrid = 'T' ; ssnd(:)%nsgn = 1. ; ssnd(:)%nct = 1
ssnd(:)%laction = .FALSE.
ssnd(:)%clgrid = 'T'
ssnd(:)%nsgn = 1.
ssnd(:)%nct = 1
ssnd(:)%dimensions = 2
! ! ------------------------- !
! ! Surface temperature !
! ! ------------------------- !
......@@ -978,6 +1094,25 @@ CONTAINS
! ! ------------------------- !
ssnd(jps_co2)%clname = 'O_CO2FLX' ; IF( TRIM(sn_snd_co2%cldes) == 'coupled' ) ssnd(jps_co2 )%laction = .TRUE.
!
#if defined key_medusa
! ! ------------------------- !
! ! MEDUSA output fields !
! ! ------------------------- !
! Surface dimethyl sulphide from Medusa
ssnd(jps_bio_dms)%clname = 'OBioDMS'
IF( TRIM(sn_snd_bio_dms%cldes) == 'medusa' ) ssnd(jps_bio_dms )%laction = .TRUE.
! Surface CO2 flux from Medusa
ssnd(jps_bio_co2)%clname = 'OBioCO2'
IF( TRIM(sn_snd_bio_co2%cldes) == 'medusa' ) ssnd(jps_bio_co2 )%laction = .TRUE.
! Surface chlorophyll from Medusa
ssnd(jps_bio_chloro)%clname = 'OBioChlo'
IF( TRIM(sn_snd_bio_chloro%cldes) == 'medusa' ) ssnd(jps_bio_chloro )%laction = .TRUE.
#endif
! ! ------------------------- !
! ! Sea surface freezing temp !
! ! ------------------------- !
......@@ -1105,10 +1240,52 @@ CONTAINS
ENDIF
ENDIF
! Initialise 1D river outflow scheme
nn_cpl_river = 1
IF ( TRIM( sn_rcv_rnf%cldes ) == 'coupled1d' ) CALL cpl_rnf_1d_init ! Coupled runoff using 1D array
! =================================================== !
! Allocate all parts of frcv used for received fields !
! =================================================== !
DO jn = 1, jprcv
IF ( srcv(jn)%laction ) THEN
SELECT CASE( srcv(jn)%dimensions )
!
CASE( 0 ) ! Scalar field
ALLOCATE( frcv(jn)%z3(1,1,1) )
CASE( 1 ) ! 1D field
ALLOCATE( frcv(jn)%z3(nn_cpl_river,1,1) )
CASE DEFAULT ! 2D (or pseudo 3D) field.
ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
END SELECT
END IF
END DO
! Allocate taum part of frcv which is used even when not received as coupling field
IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
! Allocate w10m part of frcv which is used even when not received as coupling field
IF ( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) )
! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field
IF ( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) )
IF ( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) )
! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
IF( k_ice /= 0 ) THEN
IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) )
IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) )
END IF
!
! ================================ !
! initialisation of the coupler !
! ================================ !
! There's no point initialising the coupler if we've accumulated any errors in
! coupling field definitions or settings.
IF (nstop > 0) CALL ctl_stop( 'STOP', 'sbc_cpl_init: Errors encountered in coupled field definitions' )
CALL cpl_define(jprcv, jpsnd, nn_cplmodel)
IF(ln_usecplmask) THEN
......@@ -1122,6 +1299,26 @@ CONTAINS
ENDIF
xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
!
IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
! Crude masks to separate the Antarctic and Greenland icesheets. Obviously something
! more complicated could be done if required.
greenland_icesheet_mask = 0.0
WHERE( gphit >= 0.0 ) greenland_icesheet_mask = 1.0
antarctica_icesheet_mask = 0.0
WHERE( gphit < 0.0 ) antarctica_icesheet_mask = 1.0
IF( .not. ln_rstart ) THEN
greenland_icesheet_mass = 0.0
greenland_icesheet_mass_rate_of_change = 0.0
greenland_icesheet_timelapsed = 0.0
antarctica_icesheet_mass = 0.0
antarctica_icesheet_mass_rate_of_change = 0.0
antarctica_icesheet_timelapsed = 0.0
ENDIF
ENDIF
!
IF (ln_timing) CALL timing_stop('sbc_cpl_init')
!
END SUBROUTINE sbc_cpl_init
......@@ -1181,15 +1378,26 @@ CONTAINS
LOGICAL :: llnewtx, llnewtau ! update wind stress components and module??
INTEGER :: ji, jj, jn ! dummy loop indices
INTEGER :: isec ! number of seconds since nit000 (assuming rdt did not change since nit000)
INTEGER :: ikchoix
REAL(wp), DIMENSION(jpi,jpj) :: ztx2, zty2
REAL(wp) :: zcumulneg, zcumulpos ! temporary scalars
REAL(wp) :: zcoef ! temporary scalar
LOGICAL :: ll_wrtstp ! write diagnostics?
REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3
REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient
REAL(wp) :: zgreenland_icesheet_mass_in, zantarctica_icesheet_mass_in
REAL(wp) :: zgreenland_icesheet_mass_b, zantarctica_icesheet_mass_b
REAL(wp) :: zmask_sum, zepsilon
REAL(wp) :: zzx, zzy ! temporary variables
REAL(wp) :: r1_grau ! = 1.e0 / (grav * rho0)
REAL(wp), DIMENSION(jpi,jpj) :: ztx, zty, zmsk, zemp, zqns, zqsr, zcloud_fra
!!----------------------------------------------------------------------
!
!
IF (ln_timing) CALL timing_start('sbc_cpl_rcv')
!
ll_wrtstp = ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )
!
IF( kt == nit000 ) THEN
! cannot be done in the init phase when we use agrif as cpl_freq requires that oasis_enddef is done
ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'I_QsrOce' ) + cpl_freq( 'I_QsrMix' )
......@@ -1209,7 +1417,16 @@ CONTAINS
! ! ======================================================= !
isec = ( kt - nit000 ) * NINT( rn_Dt ) ! date of exchanges
DO jn = 1, jprcv ! received fields sent by the atmosphere
IF( srcv(jn)%laction ) CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) )
IF( srcv(jn)%laction ) THEN
IF ( srcv(jn)%dimensions <= 1 ) THEN
CALL cpl_rcv_1d( jn, isec, frcv(jn)%z3, SIZE(frcv(jn)%z3), nrcvinfo(jn) )
ELSE
CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) )
END IF
END IF
END DO
! ! ========================= !
......@@ -1238,14 +1455,42 @@ CONTAINS
!
IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN ! 2 components oriented along the local grid
! ! (geographical to local grid -> rotate the components)
CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )
IF( srcv(jpr_otx2)%laction ) THEN
CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )
IF( srcv(jpr_otx1)%clgrid == 'U' .AND. (.NOT. srcv(jpr_otx2)%laction) ) THEN
! Temporary code for HadGEM3 - will be removed eventually.
! Only applies when we have only taux on U grid and tauy on V grid
!RSRH these MUST be initialised because the halos are not explicitly set
! but they're passed to repcmo and used directly in calculations, so if
! they point at junk in memory then bad things will happen!
! (You can prove this by running with preset NaNs).
ztx(:,:)=0.0
zty(:,:)=0.0
DO_2D( 0, 0, 0, 0 )
ztx(ji,jj)=0.25*vmask(ji,jj,1) &
*(frcv(jpr_otx1)%z3(ji,jj,1)+frcv(jpr_otx1)%z3(ji-1,jj,1) &
+frcv(jpr_otx1)%z3(ji,jj+1,1)+frcv(jpr_otx1)%z3(ji-1,jj+1,1))
zty(ji,jj)=0.25*umask(ji,jj,1) &
*(frcv(jpr_oty1)%z3(ji,jj,1)+frcv(jpr_oty1)%z3(ji+1,jj,1) &
+frcv(jpr_oty1)%z3(ji,jj-1,1)+frcv(jpr_oty1)%z3(ji+1,jj-1,1))
END_2D
ikchoix = 1
CALL repcmo (frcv(jpr_otx1)%z3(:,:,1),zty,ztx,frcv(jpr_oty1)%z3(:,:,1),ztx2,zty2,ikchoix)
CALL lbc_lnk ('jpr_otx1', ztx2,'U', -1. )
CALL lbc_lnk ('jpr_oty1', zty2,'V', -1. )
frcv(jpr_otx1)%z3(:,:,1)=ztx2(:,:)
frcv(jpr_oty1)%z3(:,:,1)=zty2(:,:)
ELSE
CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty )
CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )
frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:) ! overwrite 1st component on the 1st grid
IF( srcv(jpr_otx2)%laction ) THEN
CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )
ELSE
CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty )
ENDIF
frcv(jpr_oty1)%z3(:,:,1) = zty(:,:) ! overwrite 2nd component on the 2nd grid
ENDIF
frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:) ! overwrite 1st component on the 1st grid
frcv(jpr_oty1)%z3(:,:,1) = zty(:,:) ! overwrite 2nd component on the 2nd grid
ENDIF
!
IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
......@@ -1335,6 +1580,12 @@ CONTAINS
!
ENDIF
#if defined key_medusa
IF (ln_medusa) THEN
IF( srcv(jpr_atm_pco2)%laction) PCO2a_in_cpl(:,:) = frcv(jpr_atm_pco2)%z3(:,:,1)
IF( srcv(jpr_atm_dust)%laction) Dust_in_cpl(:,:) = frcv(jpr_atm_dust)%z3(:,:,1)
ENDIF
#endif
! ! ================== !
! ! atmosph. CO2 (ppm) !
! ! ================== !
......@@ -1539,8 +1790,10 @@ CONTAINS
!
IF( srcv(jpr_icb)%laction ) zqns(:,:) = zqns(:,:) - frcv(jpr_icb)%z3(:,:,1) * rLfus ! remove heat content associated to iceberg melting
!
IF( ln_mixcpl ) THEN ; qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)
ELSE ; qns(:,:) = zqns(:,:)
IF( ln_mixcpl ) THEN
qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)
ELSE
qns(:,:) = zqns(:,:)
ENDIF
! ! solar flux over the ocean (qsr)
......@@ -1559,7 +1812,78 @@ CONTAINS
IF( srcv(jpr_fice )%laction ) fr_i(:,:) = frcv(jpr_fice )%z3(:,:,1)
!
ENDIF
zepsilon = rn_iceshelf_fluxes_tolerance
IF( srcv(jpr_grnm)%laction .AND. nn_coupled_iceshelf_fluxes == 1 ) THEN
! This is a zero dimensional, single value field.
zgreenland_icesheet_mass_in = frcv(jpr_grnm)%z3(1,1,1)
greenland_icesheet_timelapsed = greenland_icesheet_timelapsed + rdt
IF( ln_iceshelf_init_atmos .AND. kt == 1 ) THEN
! On the first timestep (of an NRUN) force the ocean to ignore the icesheet masses in the ocean restart
! and take them from the atmosphere to avoid problems with using inconsistent ocean and atmosphere restarts.
zgreenland_icesheet_mass_b = zgreenland_icesheet_mass_in
greenland_icesheet_mass = zgreenland_icesheet_mass_in
ENDIF
IF( ABS( zgreenland_icesheet_mass_in - greenland_icesheet_mass ) > zepsilon ) THEN
zgreenland_icesheet_mass_b = greenland_icesheet_mass
! Only update the mass if it has increased.
IF ( (zgreenland_icesheet_mass_in - greenland_icesheet_mass) > 0.0 ) THEN
greenland_icesheet_mass = zgreenland_icesheet_mass_in
ENDIF
IF( zgreenland_icesheet_mass_b /= 0.0 ) &
& greenland_icesheet_mass_rate_of_change = ( greenland_icesheet_mass - zgreenland_icesheet_mass_b ) / greenland_icesheet_timelapsed
greenland_icesheet_timelapsed = 0.0_wp
ENDIF
IF(lwp .AND. ll_wrtstp) THEN
WRITE(numout,*) 'Greenland icesheet mass (kg) read in is ', zgreenland_icesheet_mass_in
WRITE(numout,*) 'Greenland icesheet mass (kg) used is ', greenland_icesheet_mass
WRITE(numout,*) 'Greenland icesheet mass rate of change (kg/s) is ', greenland_icesheet_mass_rate_of_change
WRITE(numout,*) 'Greenland icesheet seconds lapsed since last change is ', greenland_icesheet_timelapsed
ENDIF
ELSE IF ( nn_coupled_iceshelf_fluxes == 2 ) THEN
greenland_icesheet_mass_rate_of_change = rn_greenland_total_fw_flux
ENDIF
! ! land ice masses : Antarctica
IF( srcv(jpr_antm)%laction .AND. nn_coupled_iceshelf_fluxes == 1 ) THEN
! This is a zero dimensional, single value field.
zantarctica_icesheet_mass_in = frcv(jpr_antm)%z3(1,1,1)
antarctica_icesheet_timelapsed = antarctica_icesheet_timelapsed + rdt
IF( ln_iceshelf_init_atmos .AND. kt == 1 ) THEN
! On the first timestep (of an NRUN) force the ocean to ignore the icesheet masses in the ocean restart
! and take them from the atmosphere to avoid problems with using inconsistent ocean and atmosphere restarts.
zantarctica_icesheet_mass_b = zantarctica_icesheet_mass_in
antarctica_icesheet_mass = zantarctica_icesheet_mass_in
ENDIF
IF( ABS( zantarctica_icesheet_mass_in - antarctica_icesheet_mass ) > zepsilon ) THEN
zantarctica_icesheet_mass_b = antarctica_icesheet_mass
! Only update the mass if it has increased.
IF ( (zantarctica_icesheet_mass_in - antarctica_icesheet_mass) > 0.0 ) THEN
antarctica_icesheet_mass = zantarctica_icesheet_mass_in
END IF
IF( zantarctica_icesheet_mass_b /= 0.0 ) &
& antarctica_icesheet_mass_rate_of_change = ( antarctica_icesheet_mass - zantarctica_icesheet_mass_b ) / antarctica_icesheet_timelapsed
antarctica_icesheet_timelapsed = 0.0_wp
ENDIF
IF(lwp .AND. ll_wrtstp) THEN
WRITE(numout,*) 'Antarctica icesheet mass (kg) read in is ', zantarctica_icesheet_mass_in
WRITE(numout,*) 'Antarctica icesheet mass (kg) used is ', antarctica_icesheet_mass
WRITE(numout,*) 'Antarctica icesheet mass rate of change (kg/s) is ', antarctica_icesheet_mass_rate_of_change
WRITE(numout,*) 'Antarctica icesheet seconds lapsed since last change is ', antarctica_icesheet_timelapsed
ENDIF
ELSE IF ( nn_coupled_iceshelf_fluxes == 2 ) THEN
antarctica_icesheet_mass_rate_of_change = rn_antarctica_total_fw_flux
ENDIF
!
IF (ln_timing) CALL timing_stop('sbc_cpl_rcv')
END SUBROUTINE sbc_cpl_rcv
......@@ -1759,17 +2083,32 @@ CONTAINS
!!----------------------------------------------------------------------
!
#if defined key_si3 || defined key_cice
!
IF( kt == nit000 ) THEN
! allocate ice fractions from last coupling time here and not in sbc_cpl_init because of jpl
IF( .NOT.ALLOCATED(a_i_last_couple) ) ALLOCATE( a_i_last_couple(jpi,jpj,jpl) )
! initialize to a_i for the 1st time step
a_i_last_couple(:,:,:) = a_i(:,:,:)
ENDIF
!
IF( ln_mixcpl ) zmsk(:,:) = 1. - xcplmask(:,:,0)
ziceld(:,:) = 1._wp - picefr(:,:)
zcptn (:,:) = rcp * sst_m(:,:)
! ! ========================= !
SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) ) ! ice topmelt and botmelt !
! ! ========================= !
CASE ('coupled')
IF (ln_scale_ice_flux) THEN
WHERE( a_i(:,:,:) > 1.e-10_wp )
qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
ELSEWHERE
qml_ice(:,:,:) = 0.0_wp
qcn_ice(:,:,:) = 0.0_wp
END WHERE
ELSE
qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
ENDIF
END SELECT
!
! ! ========================= !
! ! freshwater budget ! (emp_tot)
......@@ -1783,7 +2122,9 @@ CONTAINS
CASE( 'conservative' ) ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
zsprecip(:,:) = frcv(jpr_snow)%z3(:,:,1) ! May need to ensure positive here
ztprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:) ! May need to ensure positive here
zemp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
IF (.not. ln_couple_ocean_evap ) THEN
zemp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
END IF
CASE( 'oce and ice' ) ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
zemp_tot(:,:) = ziceld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + picefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) * picefr(:,:)
......@@ -1791,7 +2132,9 @@ CONTAINS
ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)
CASE( 'none' ) ! Not available as for now: needs additional coding below when computing zevap_oce
! ! since fields received are not defined with none option
CALL ctl_stop('STOP', 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_emp value in namelist namsbc_cpl')
CALL ctl_stop( 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_emp value in namelist namsbc_cpl' )
CASE default ! Default
CALL ctl_stop( 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_emp value in namelist namsbc_cpl' )
END SELECT
! --- evaporation over ice (kg/m2/s) --- !
......@@ -1837,19 +2180,30 @@ CONTAINS
! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- !
zemp_ice(:,:) = zemp_ice(:,:) + zsprecip(:,:) * ( picefr(:,:) - zsnw(:,:) ) ! emp_ice = A * sublimation - zsnw * sprecip
zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:) ! emp_oce = emp_tot - emp_ice
! --- evaporation over ocean (used later for qemp) --- !
zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:)
IF ( ln_couple_ocean_evap ) THEN
zemp_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_rain)%z3(:,:,1) & !Ocean evap minus rain (as all rain goes straight to ocean in GC5)
& - zsprecip(:,:) * ( 1._wp - zsnw(:,:) ) !subtract snow in leads after correction for blowing snow
zemp_tot(:,:) = zemp_oce(:,:) + zemp_ice(:,:)
zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1)
ELSE
zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:) ! emp_oce = emp_tot - emp_ice
! --- evaporation over ocean (used later for qemp) --- !
zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:)
END IF
! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0
! therefore, sublimation is not redistributed over the ice categories when no subgrid scale fluxes are provided by atm.
zdevap_ice(:,:) = 0._wp
! --- Continental fluxes --- !
IF( srcv(jpr_rnf)%laction ) THEN ! runoffs (included in emp later on)
IF( srcv(jpr_rnf)%laction ) THEN ! 2D runoffs (included in emp later on)
rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
ENDIF
IF( srcv(jpr_rnf_1d)%laction ) THEN ! 1D runoff
CALL cpl_rnf_1d_to_2d(frcv(jpr_rnf_1d)%z3(:,:,:))
ENDIF
IF( srcv(jpr_cal)%laction ) THEN ! calving (put in emp_tot and emp_oce)
zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
zemp_oce(:,:) = zemp_oce(:,:) - frcv(jpr_cal)%z3(:,:,1)
......@@ -1929,23 +2283,7 @@ CONTAINS
!! IF( srcv(jpr_rnf)%laction ) CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1) ) ! runoff
!! IF( srcv(jpr_isf)%laction ) CALL iom_put( 'iceshelf_cea', fwfisf(:,:) * tmask(:,:,1) ) ! iceshelf
!
! ! ========================= !
SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) ) ! ice topmelt and botmelt !
! ! ========================= !
CASE ('coupled')
IF (ln_scale_ice_flux) THEN
WHERE( a_i(:,:,:) > 1.e-10_wp )
qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
ELSEWHERE
qml_ice(:,:,:) = 0.0_wp
qcn_ice(:,:,:) = 0.0_wp
END WHERE
ELSE
qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
ENDIF
END SELECT
!
! ! ========================= !
SELECT CASE( TRIM( sn_rcv_qns%cldes ) ) ! non solar heat fluxes ! (qns)
......@@ -2009,8 +2347,10 @@ CONTAINS
! --- non solar flux over ocean --- !
! note: ziceld cannot be = 0 since we limit the ice concentration to amax
zqns_oce = 0._wp
WHERE( ziceld /= 0._wp ) zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / ziceld(:,:)
WHERE( ziceld /= 0._wp ) zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / ziceld(:,:)
! Heat content per unit mass of snow (J/kg)
WHERE( SUM( a_i, dim=3 ) > 1.e-10 ) ; zcptsnw(:,:) = rcpi * SUM( (tn_ice - rt0) * a_i, dim=3 ) / SUM( a_i, dim=3 )
ELSEWHERE ; zcptsnw(:,:) = zcptn(:,:)
......@@ -2023,7 +2363,12 @@ CONTAINS
! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
DO jl = 1, jpl
zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * rcpi ) but atm. does not take it into account
! zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * rcpi ) but atm. does not take it into account
! How much enthalpy is stored in sublimating snow and ice. At this stage we don't know if it is snow or ice that is
! sublimating so we will use the combined snow and ice layer temperature t1_ice.
zqevap_ice(:,:,jl) = -zevap_ice(:,:,jl) * ( ( rt0 - t1_ice(:,:,jl) ) * rcpi + rLfus )
END DO
! --- heat flux associated with emp (W/m2) --- !
......@@ -2041,6 +2386,7 @@ CONTAINS
IF( ln_mixcpl ) THEN
qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) + zqns_tot(:,:)* zmsk(:,:)
qns_oce(:,:) = qns_oce(:,:) * xcplmask(:,:,0) + zqns_oce(:,:)* zmsk(:,:)
DO jl=1,jpl
qns_ice (:,:,jl) = qns_ice (:,:,jl) * xcplmask(:,:,0) + zqns_ice (:,:,jl)* zmsk(:,:)
qevap_ice(:,:,jl) = qevap_ice(:,:,jl) * xcplmask(:,:,0) + zqevap_ice(:,:,jl)* zmsk(:,:)
......@@ -2085,9 +2431,14 @@ CONTAINS
IF ( srcv(jpr_icb)%laction ) CALL iom_put('hflx_icb_cea' , - frcv(jpr_icb)%z3(:,:,1) * rLfus ) ! latent heat from icebergs melting
IF ( iom_use('hflx_rain_cea') ) & ! heat flux from rain (cell average)
& CALL iom_put('hflx_rain_cea' , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )
IF ( iom_use('hflx_evap_cea') ) & ! heat flux from evap (cell average)
& CALL iom_put('hflx_evap_cea' , ( frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:) ) &
& * zcptn(:,:) * tmask(:,:,1) )
IF ( ln_couple_ocean_evap ) THEN
IF ( iom_use( 'hflx_evap_cea') ) CALL iom_put('hflx_evap_cea' , frcv(jpr_tevp)%z3(:,:,1) * zcptn(:,:) * tmask(:,:,1) ) ! heat flux from evap (cell average)
ELSE
IF ( iom_use( 'hflx_evap_cea') ) CALL iom_put('hflx_evap_cea' , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) ) &
& * zcptn(:,:) * tmask(:,:,1) ) ! heat flux from evap (cell average)
END IF
IF ( iom_use('hflx_prec_cea') ) & ! heat flux from all precip (cell avg)
& CALL iom_put('hflx_prec_cea' , sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) &
& + ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )
......@@ -2214,29 +2565,30 @@ CONTAINS
!
ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN !== conduction flux as surface forcing ==!
!
!! SELECT CASE( TRIM( sn_rcv_qtrice%cldes ) )
!! !
!! ! ! ===> here we receive the qtr_ice_top array from the coupler
!! CASE ('coupled')
!! IF (ln_scale_ice_flux) THEN
!! WHERE( a_i(:,:,:) > 1.e-10_wp )
!! zqtr_ice_top(:,:,:) = frcv(jpr_qtrice)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
!! ELSEWHERE
!! zqtr_ice_top(:,:,:) = 0.0_wp
!! ENDWHERE
!! ELSE
!! zqtr_ice_top(:,:,:) = frcv(jpr_qtrice)%z3(:,:,:)
!! ENDIF
!!
!! ! Add retrieved transmitted solar radiation onto the ice and total solar radiation
!! zqsr_ice(:,:,:) = zqsr_ice(:,:,:) + zqtr_ice_top(:,:,:)
!! zqsr_tot(:,:) = zqsr_tot(:,:) + SUM( zqtr_ice_top(:,:,:) * a_i(:,:,:), dim=3 )
!!
!! ! if we are not getting this data from the coupler then assume zero (fully opaque ice)
!! CASE ('none')
zqtr_ice_top(:,:,:) = 0._wp
!! END SELECT
!
SELECT CASE( TRIM( sn_rcv_qtr%cldes ) )
!
! ! ===> here we receive the qtr_ice_top array from the coupler
CASE ('coupled')
IF (ln_scale_ice_flux) THEN
WHERE( a_i(:,:,:) > 0.0_wp ) zqtr_ice_top(:,:,:) = MAX(0._wp, frcv(jpr_qtr)%z3(:,:,:)) * a_i_last_couple(:,:,:) / a_i(:,:,:)
WHERE( a_i(:,:,:) <= 0.0_wp ) zqtr_ice_top(:,:,:) = 0.0_wp
ELSE
zqtr_ice_top(:,:,:) = MAX(0._wp, frcv(jpr_qtr)%z3(:,:,:))
ENDIF
! Add retrieved transmitted solar radiation onto the ice and total solar radiation
zqsr_ice(:,:,:) = zqsr_ice(:,:,:) + zqtr_ice_top(:,:,:)
zqsr_tot(:,:) = zqsr_tot(:,:) + SUM( zqtr_ice_top(:,:,:) * a_i(:,:,:), dim=3 )
! if we are not getting this data from the coupler then assume zero (fully opaque ice)
CASE ('none')
zqtr_ice_top(:,:,:) = 0._wp
CASE default
CALL ctl_stop( 'sbc_cpl_ice_flx: Invalid value for sn_rcv_qtr%cldes' )
END SELECT
!
ENDIF
IF( ln_mixcpl ) THEN
......@@ -2260,6 +2612,44 @@ CONTAINS
IF( ln_mixcpl ) THEN ; qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) + zqsr_oce(:,:)* zmsk(:,:)
ELSE ; qsr_oce(:,:) = zqsr_oce(:,:) ; ENDIF
IF( ln_mixcpl ) THEN
qsr_tot(:,:) = qsr(:,:) * ziceld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 ) ! total flux from blk
qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) + zqsr_tot(:,:)* zmsk(:,:)
DO jl = 1, jpl
qsr_ice(:,:,jl) = qsr_ice(:,:,jl) * xcplmask(:,:,0) + zqsr_ice(:,:,jl)* zmsk(:,:)
END DO
ELSE
qsr_tot(:,: ) = zqsr_tot(:,: )
qsr_ice(:,:,:) = zqsr_ice(:,:,:)
ENDIF
! ! ========================= !
SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) ) ! d(qns)/dt !
! ! ========================= !
CASE ('coupled')
IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
zdqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
ELSE
! Set all category values equal for the moment
DO jl=1,jpl
zdqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
ENDDO
ENDIF
CASE( 'none' )
zdqns_ice(:,:,:) = 0._wp
CASE default
CALL ctl_stop( 'sbc_cpl_ice_flx: Invalid value for sn_rcv_dqnsdt%cldes' )
END SELECT
IF( ln_mixcpl ) THEN
DO jl=1,jpl
dqns_ice(:,:,jl) = dqns_ice(:,:,jl) * xcplmask(:,:,0) + zdqns_ice(:,:,jl) * zmsk(:,:)
ENDDO
ELSE
dqns_ice(:,:,:) = zdqns_ice(:,:,:)
ENDIF
! ! ================== !
! ! ice skin temp. !
! ! ================== !
......@@ -2297,13 +2687,18 @@ CONTAINS
INTEGER, INTENT(in) :: kt
INTEGER, INTENT(in) :: Kbb, Kmm ! ocean model time level index
!
INTEGER :: ji, jj, jl ! dummy loop indices
INTEGER :: ji, jj, jl ! dummy loop indices
INTEGER :: ikchoix
INTEGER :: isec, info ! local integer
REAL(wp) :: zumax, zvmax
REAL(wp), DIMENSION(jpi,jpj) :: zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztmp3, ztmp4
REAL(wp), DIMENSION(jpi,jpj) :: ztmp5, ztmp6 ! RSRH temporary work arrays
! to avoid intent conflicts in repcmo calls
!!----------------------------------------------------------------------
!
IF (ln_timing) CALL timing_start('sbc_cpl_snd')
isec = ( kt - nit000 ) * NINT( rn_Dt ) ! date of exchanges
info = OASIS_idle
......@@ -2317,12 +2712,15 @@ CONTAINS
ztmp1(:,:) = ts(:,:,1,jp_tem,Kmm) ! send temperature as it is (potential or conservative) -> use of l_useCT on the received part
ELSE
! we must send the surface potential temperature
IF( l_useCT ) THEN ; ztmp1(:,:) =eos_pt_from_ct( CASTSP(ts(:,:,1,jp_tem,Kmm)), CASTSP(ts(:,:,1,jp_sal,Kmm)) )
ELSE ; ztmp1(:,:) = ts(:,:,1,jp_tem,Kmm)
IF( l_useCT ) THEN
ztmp1(:,:) =eos_pt_from_ct( CASTSP(ts(:,:,1,jp_tem,Kmm)), CASTSP(ts(:,:,1,jp_sal,Kmm)) )
ELSE
ztmp1(:,:) = ts(:,:,1,jp_tem,Kmm)
ENDIF
!
SELECT CASE( sn_snd_temp%cldes)
CASE( 'oce only' ) ; ztmp1(:,:) = ztmp1(:,:) + rt0
CASE( 'oce only' )
ztmp1(:,:) = ztmp1(:,:) + rt0
CASE( 'oce and ice' ) ; ztmp1(:,:) = ztmp1(:,:) + rt0
SELECT CASE( sn_snd_temp%clcat )
CASE( 'yes' )
......@@ -2564,18 +2962,37 @@ CONTAINS
! ! ------------------------- !
! ! CO2 flux from PISCES !
! ! ------------------------- !
IF( ssnd(jps_co2)%laction .AND. l_co2cpl ) THEN
ztmp1(:,:) = oce_co2(:,:) * 1000. ! conversion in molC/m2/s
CALL cpl_snd( jps_co2, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ) , info )
IF( ssnd(jps_co2)%laction .AND. l_co2cpl ) THEN
ztmp1(:,:) = oce_co2(:,:) * 1000. ! conversion in molC/m2/s
CALL cpl_snd( jps_co2, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ) , info )
ENDIF
#if defined key_medusa
!
IF (ln_medusa) THEN
! ! ---------------------------------------------- !
! ! CO2 flux, DMS and chlorophyll from MEDUSA !
! ! ---------------------------------------------- !
IF ( ssnd(jps_bio_co2)%laction ) THEN
CALL cpl_snd( jps_bio_co2, isec, RESHAPE( CO2Flux_out_cpl, (/jpi,jpj,1/) ), info )
ENDIF
IF ( ssnd(jps_bio_dms)%laction ) THEN
CALL cpl_snd( jps_bio_dms, isec, RESHAPE( DMS_out_cpl, (/jpi,jpj,1/) ), info )
ENDIF
IF ( ssnd(jps_bio_chloro)%laction ) THEN
CALL cpl_snd( jps_bio_chloro, isec, RESHAPE( chloro_out_cpl, (/jpi,jpj,1/) ), info )
ENDIF
ENDIF
!
#endif
! ! ------------------------- !
IF( ssnd(jps_ocx1)%laction ) THEN ! Surface current !
! ! ------------------------- !
!
! j+1 j -----V---F
! surface velocity always sent from T point ! |
! j | T U
! [except for HadGEM3] j | T U
! | |
! j j-1 -I-------|
! (for I) | |
......@@ -2587,11 +3004,19 @@ CONTAINS
ELSE
SELECT CASE( TRIM( sn_snd_crt%cldes ) )
CASE( 'oce only' ) ! C-grid ==> T
DO_2D( 0, 0, 0, 0 )
zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Kmm) + uu(ji-1,jj ,1,Kmm) )
zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Kmm) + vv(ji ,jj-1,1,Kmm) )
END_2D
CASE( 'weighted oce and ice' ) ! Ocean and Ice on C-grid ==> T
IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN
DO_2D( 0, 0, 0, 0 )
zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Kmm) + uu(ji-1,jj ,1,Kmm) )
zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Kmm) + vv(ji ,jj-1,1,Kmm) )
END_2D
ELSE
! Temporarily Changed for UKV
DO_2D( 0, 0, 0, 0 )
zotx1(ji,jj) = uu(ji,jj,1,Kmm)
zoty1(ji,jj) = vv(ji,jj,1,Kmm)
END_2D
ENDIF
CASE( 'weighted oce and ice' ) ! Ocean and Ice on C-grid ==> T
DO_2D( 0, 0, 0, 0 )
zotx1(ji,jj) = 0.5 * ( uu (ji,jj,1,Kmm) + uu (ji-1,jj ,1,Kmm) ) * zfr_l(ji,jj)
zoty1(ji,jj) = 0.5 * ( vv (ji,jj,1,Kmm) + vv (ji ,jj-1,1,Kmm) ) * zfr_l(ji,jj)
......@@ -2607,6 +3032,7 @@ CONTAINS
& + 0.5 * ( v_ice(ji,jj ) + v_ice(ji ,jj-1 ) ) * fr_i(ji,jj)
END_2D
END SELECT
CALL lbc_lnk( 'sbccpl', zotx1, ssnd(jps_ocx1)%clgrid, -1.0_wp, zoty1, ssnd(jps_ocy1)%clgrid, -1.0_wp )
!
ENDIF
......@@ -2614,15 +3040,44 @@ CONTAINS
!
IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN ! Rotation of the components
! ! Ocean component
CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 ) ! 1st component
CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 ) ! 2nd component
zotx1(:,:) = ztmp1(:,:) ! overwrite the components
zoty1(:,:) = ztmp2(:,:)
IF( ssnd(jps_ivx1)%laction ) THEN ! Ice component
CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 ) ! 1st component
CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 ) ! 2nd component
zitx1(:,:) = ztmp1(:,:) ! overwrite the components
zity1(:,:) = ztmp2(:,:)
IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN
CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 ) ! 1st component
CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 ) ! 2nd component
zotx1(:,:) = ztmp1(:,:) ! overwrite the components
zoty1(:,:) = ztmp2(:,:)
IF( ssnd(jps_ivx1)%laction ) THEN ! Ice component
CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 ) ! 1st component
CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 ) ! 2nd component
zitx1(:,:) = ztmp1(:,:) ! overwrite the components
zity1(:,:) = ztmp2(:,:)
ENDIF
ELSE
! Temporary code for HadGEM3 - will be removed eventually.
! Only applies when we want uvel on U grid and vvel on V grid
! Rotate U and V onto geographic grid before sending.
DO_2D( 0, 0, 0, 0 )
ztmp1(ji,jj)=0.25*vmask(ji,jj,1) &
*(zotx1(ji,jj)+zotx1(ji-1,jj) &
+zotx1(ji,jj+1)+zotx1(ji-1,jj+1))
ztmp2(ji,jj)=0.25*umask(ji,jj,1) &
*(zoty1(ji,jj)+zoty1(ji+1,jj) &
+zoty1(ji,jj-1)+zoty1(ji+1,jj-1))
END_2D
ikchoix = -1
! zotx1 and zoty1 are input only to repcmo while ztmp5 and ztmp6
! are the newly calculated (output) values.
! Don't make the mistake of using zotx1 and zoty1 twice in this
! call for both input and output fields since it creates INTENT
! conflicts.
CALL repcmo (zotx1,ztmp2,ztmp1,zoty1,ztmp5,ztmp6,ikchoix)
zotx1(:,:)=ztmp5(:,:)
zoty1(:,:)=ztmp6(:,:)
! Ensure any N fold and wrap columns are updated.
CALL lbc_lnk( 'sbccpl', zotx1, ssnd(jps_ocx1)%clgrid, -1.0_wp, zoty1, ssnd(jps_ocy1)%clgrid, -1.0_wp )
ENDIF
ENDIF
!
......@@ -2783,6 +3238,7 @@ CONTAINS
IF( ssnd(jps_sstfrz)%laction ) CALL cpl_snd( jps_sstfrz, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info)
#endif
!
IF (ln_timing) CALL timing_stop('sbc_cpl_snd')
END SUBROUTINE sbc_cpl_snd
!!======================================================================
......
src/OCE/SBC/sbcrnf.F90
View file @ 672ff863
......@@ -20,6 +20,7 @@ MODULE sbcrnf
USE sbc_oce ! surface boundary condition variables
USE eosbn2 ! Equation Of State
USE closea, ONLY: l_clo_rnf, clo_rnf ! closed seas
USE isf_oce
!
USE in_out_manager ! I/O manager
USE fldread ! read input field at current time step
......@@ -40,8 +41,10 @@ MODULE sbcrnf
LOGICAL :: ln_rnf_depth_ini !: depth river runoffs computed at the initialisation
REAL(wp) :: rn_rnf_max !: maximum value of the runoff climatologie (ln_rnf_depth_ini =T)
REAL(wp) :: rn_dep_max !: depth over which runoffs is spread (ln_rnf_depth_ini =T)
REAL(wp) :: tot_flux !: total iceberg flux (temporary variable)
INTEGER :: nn_rnf_depth_file !: create (=1) a runoff depth file or not (=0)
LOGICAL , PUBLIC :: ln_rnf_icb !: iceberg flux is specified in a file
LOGICAL :: ln_icb_mass !: Scale iceberg flux to match FW flux from coupled model
LOGICAL :: ln_rnf_tem !: temperature river runoffs attribute specified in a file
LOGICAL , PUBLIC :: ln_rnf_sal !: salinity river runoffs attribute specified in a file
TYPE(FLD_N) , PUBLIC :: sn_rnf !: information about the runoff file to be read
......@@ -119,8 +122,8 @@ CONTAINS
!
IF( .NOT. l_rnfcpl ) THEN
CALL fld_read ( kt, nn_fsbc, sf_rnf ) ! Read Runoffs data and provide it at kt ( runoffs + iceberg )
IF( ln_rnf_icb ) CALL fld_read ( kt, nn_fsbc, sf_i_rnf ) ! idem for iceberg flux if required
ENDIF
IF( ln_rnf_icb ) CALL fld_read ( kt, nn_fsbc, sf_i_rnf ) ! idem for iceberg flux if required
IF( ln_rnf_tem ) CALL fld_read ( kt, nn_fsbc, sf_t_rnf ) ! idem for runoffs temperature if required
IF( ln_rnf_sal ) CALL fld_read ( kt, nn_fsbc, sf_s_rnf ) ! idem for runoffs salinity if required
!
......@@ -138,14 +141,57 @@ CONTAINS
CALL iom_put( 'hflx_icb_cea' , -fwficb(:,:) * rLfus ) ! output Heat Flux into Sea Water due to Iceberg Thermodynamics -->
ENDIF
ENDIF
!
IF( ln_rnf_icb ) THEN
fwficb(:,:) = rn_rfact * ( sf_i_rnf(1)%fnow(:,:,1) ) * tmask(:,:,1) ! updated runoff value at time step kt
IF( ln_icb_mass ) THEN
! Modify the Iceberg FW flux to be consistent with the change in
! mass of the Antarctic/Greenland ice sheet for FW conservation in
! coupled model. This isn't perfect as FW flux will go into ocean at
! wrong time of year but more important to maintain FW balance
tot_flux = SUM(fwficb(:,:)*e1e2t(:,:)*tmask_i(:,:)*greenland_icesheet_mask(:,:)) !Need to multiply by area to convert to kg/s
IF( lk_mpp ) CALL mpp_sum( 'icbclv', tot_flux )
IF( tot_flux > rsmall ) THEN
WHERE( greenland_icesheet_mask(:,:) == 1.0 ) &
& fwficb(:,:) = fwficb(:,:) * greenland_icesheet_mass_rate_of_change * rn_greenland_calving_fraction &
& / tot_flux
ELSE IF( rn_greenland_calving_fraction < rsmall ) THEN
WHERE( greenland_icesheet_mask(:,:) == 1.0 ) fwficb(:,:) = 0.0
ELSE
CALL CTL_STOP('STOP', 'No iceberg runoff data read in for Greenland. Check input file or set rn_greenland_calving_fraction=0.0')
ENDIF
tot_flux = SUM(fwficb(:,:)*e1e2t(:,:)*tmask_i(:,:)*antarctica_icesheet_mask(:,:))
IF( lk_mpp ) CALL mpp_sum( 'icbclv', tot_flux )
IF( tot_flux > rsmall ) THEN
WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
& fwficb(:,:) = fwficb(:,:) * antarctica_icesheet_mass_rate_of_change * rn_antarctica_calving_fraction &
& / (tot_flux + 1.0e-10_wp )
ELSE IF( rn_antarctica_calving_fraction < rsmall ) THEN
WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) fwficb(:,:) = 0.0
ELSE
CALL CTL_STOP('STOP', 'No iceberg runoff data read in for Antarctica. Check input file or set rn_antarctica_calving_fraction=0.0')
ENDIF
ENDIF
CALL lbc_lnk('sbcrnf',rnf,'T',1._wp,fwficb,'T',1._wp)
CALL iom_put( 'berg_melt' , fwficb(:,:) ) ! output iceberg flux
CALL iom_put( 'hflx_icb_cea' , fwficb(:,:) * rLfus ) ! output Heat Flux into Sea Water due to Iceberg Thermodynamics -->
rnf(:,:) = rnf(:,:) + fwficb(:,:) ! fwficb isn't used anywhere else so add it to runoff here
qns_tot(:,:) = qns_tot(:,:) - fwficb(:,:) * rLfus ! TG: I think this is correct
ENDIF !
! ! set temperature & salinity content of runoffs
IF( ln_rnf_tem ) THEN ! use runoffs temperature data
rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rho0
CALL eos_fzp( sss_m(:,:), ztfrz(:,:) )
WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp ) ! if missing data value use SST as runoffs temperature
rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rho0
rnf_tsc(:,:,jp_tem) = MAX( sst_m(:,:), 0.0_wp ) * rnf(:,:) * r1_rho0
END WHERE
WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )
! where fwf comes from melting of ice shelves or iceberg
rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rho0 - rnf(:,:) * rLfusisf * r1_rho0_rcp
END WHERE
ELSE ! use SST as runoffs temperature
!CEOD River is fresh water so must at least be 0 unless we consider ice
rnf_tsc(:,:,jp_tem) = MAX( sst_m(:,:), 0.0_wp ) * rnf(:,:) * r1_rho0
......@@ -158,11 +204,15 @@ CONTAINS
IF( iom_use('sflx_rnf_cea') ) CALL iom_put( 'sflx_rnf_cea', rnf_tsc(:,:,jp_sal) * rho0 ) ! output runoff salt flux (g/m2/s)
ENDIF
!
! Make sure rnf is up to date!
call lbc_lnk('rnf_div', rnf, 'T', 1.0)
! ! ---------------------------------------- !
IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 !
! ! ---------------------------------------- !
IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN !* Restart: read in restart file
IF(lwp) WRITE(numout,*) ' nit000-1 runoff forcing fields red in the restart file', lrxios
IF(lwp) WRITE(numout,*) ' nit000-1 runoff forcing fields read in the restart file', lrxios
CALL iom_get( numror, jpdom_auto, 'rnf_b' , rnf_b ) ! before runoff
CALL iom_get( numror, jpdom_auto, 'rnf_hc_b', rnf_tsc_b(:,:,jp_tem) ) ! before heat content of runoff
CALL iom_get( numror, jpdom_auto, 'rnf_sc_b', rnf_tsc_b(:,:,jp_sal) ) ! before salinity content of runoff
......@@ -259,7 +309,7 @@ CONTAINS
REAL(wp) :: zacoef
REAL(wp), DIMENSION(jpi,jpj,2) :: zrnfcl
!!
NAMELIST/namsbc_rnf/ cn_dir , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal, ln_rnf_icb, &
NAMELIST/namsbc_rnf/ cn_dir , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal, ln_rnf_icb, ln_icb_mass, &
& sn_rnf, sn_cnf , sn_i_rnf, sn_s_rnf , sn_t_rnf , sn_dep_rnf, &
& ln_rnf_mouth , rn_hrnf , rn_avt_rnf, rn_rfact, &
& ln_rnf_depth_ini , rn_dep_max , rn_rnf_max, nn_rnf_depth_file
......@@ -318,21 +368,22 @@ CONTAINS
IF( sn_rnf%ln_tint ) ALLOCATE( sf_rnf(1)%fdta(jpi,jpj,1,2) )
CALL fld_fill( sf_rnf, (/ sn_rnf /), cn_dir, 'sbc_rnf_init', 'read runoffs data', 'namsbc_rnf', no_print )
!
IF( ln_rnf_icb ) THEN ! Create (if required) sf_i_rnf structure
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' iceberg flux read in a file'
ALLOCATE( sf_i_rnf(1), STAT=ierror )
IF( ierror > 0 ) THEN
CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_i_rnf structure' ) ; RETURN
ENDIF
ALLOCATE( sf_i_rnf(1)%fnow(jpi,jpj,1) )
IF( sn_i_rnf%ln_tint ) ALLOCATE( sf_i_rnf(1)%fdta(jpi,jpj,1,2) )
CALL fld_fill (sf_i_rnf, (/ sn_i_rnf /), cn_dir, 'sbc_rnf_init', 'read iceberg flux data', 'namsbc_rnf' )
ELSE
fwficb(:,:) = 0._wp
ENDIF
ENDIF
IF( ln_rnf_icb ) THEN ! Create (if required) sf_i_rnf structure
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' iceberg flux read in a file'
ALLOCATE( sf_i_rnf(1), STAT=ierror )
IF( ierror > 0 ) THEN
CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_i_rnf structure' ) ; RETURN
ENDIF
ALLOCATE( sf_i_rnf(1)%fnow(jpi,jpj,1) )
IF( sn_i_rnf%ln_tint ) ALLOCATE( sf_i_rnf(1)%fdta(jpi,jpj,1,2) )
CALL fld_fill (sf_i_rnf, (/ sn_i_rnf /), cn_dir, 'sbc_rnf_init', 'read iceberg flux data', 'namsbc_rnf' )
ELSE
fwficb(:,:) = 0._wp
ENDIF
!
IF( ln_rnf_tem ) THEN ! Create (if required) sf_t_rnf structure
IF(lwp) WRITE(numout,*)
......
src/OCE/SBC/sbcssm.F90
View file @ 672ff863
......@@ -59,13 +59,28 @@ CONTAINS
INTEGER :: ji, jj ! loop index
REAL(wp) :: zcoef, zf_sbc ! local scalar
REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts
!!---------------------------------------------------------------------
CHARACTER(len=4),SAVE :: stype
!!---------------------------------------------------------------------
IF( kt == nit000 ) THEN
IF( ln_TEOS10 ) THEN
stype='abs' ! teos-10: using absolute salinity (sst is converted to potential temperature for the surface module)
ELSE IF( ln_SEOS ) THEN
stype='seos' ! seos: using Simplified Equation of state (sst is converted to potential temperature for the surface module)
ELSE
stype='pra' ! eos-80: using practical salinity
ENDIF
ENDIF
!
! !* surface T-, U-, V- ocean level variables (T, S, depth, velocity)
zts(:,:,jp_tem) = ts(:,:,1,jp_tem,Kmm)
zts(:,:,jp_sal) = ts(:,:,1,jp_sal,Kmm)
!
! ! ---------------------------------------- !
! !===>>> CAUTION: lbc_lnk is required on fraqsr_lev since sea ice computes on the full domain
! ! otherwise restartability and reproducibility are broken
! ! computed in traqsr only on the inner domain
CALL lbc_lnk( 'sbc_ssm', fraqsr_1lev(:,:), 'T', 1._wp )
!
! ! ---------------------------------------- !
IF( nn_fsbc == 1 ) THEN ! Instantaneous surface fields !
! ! ---------------------------------------- !
ssu_m(:,:) = uu(:,:,1,Kbb)
......@@ -171,8 +186,8 @@ CONTAINS
IF( MOD( kt - 1 , nn_fsbc ) == 0 ) THEN ! Mean value at each nn_fsbc time-step !
CALL iom_put( 'ssu_m', ssu_m )
CALL iom_put( 'ssv_m', ssv_m )
CALL iom_put( 'sst_m', sst_m )
CALL iom_put( 'sss_m', sss_m )
CALL iom_put( 'sst_m_pot', sst_m )
CALL iom_put( 'sss_m_'//stype, sss_m )
CALL iom_put( 'ssh_m', ssh_m )
CALL iom_put( 'e3t_m', e3t_m )
CALL iom_put( 'frq_m', frq_m )
......
src/OCE/TRA/eosbn2.F90
View file @ 672ff863
......@@ -45,6 +45,7 @@ MODULE eosbn2
!
USE in_out_manager ! I/O manager
USE lib_mpp ! for ctl_stop
USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
USE prtctl ! Print control
USE timing ! Timing
......
src/OCE/TRA/traadv_fct.F90
View file @ 672ff863
......@@ -87,6 +87,7 @@ CONTAINS
REAL(wp) :: zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u ! - -
REAL(wp) :: zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v ! - -
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwi, zwz, ztu, ztv, zltu, zltv, ztw
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwi_in, ztw_in !temp arrays to avoid intent conflicts
REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zwx, zwy
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz, zptry
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zwinf, zwdia, zwsup
......@@ -194,7 +195,9 @@ CONTAINS
END_3D
IF ( ll_zAimp ) THEN
CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 )
! We need to use separate copies of zwi to avoid intent conflicts!
zwi_in(:,:,:) = zwi(:,:,:)
CALL tridia_solver( zwdia, zwsup, zwinf, zwi_in, zwi , 0 )
!
ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Interior value ( multiplied by wmask)
......@@ -314,7 +317,10 @@ CONTAINS
ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
END_3D
!
CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
! We need to use separate copies of ztw to avoid intent conflicts!
ztw_in(:,:,:) = ztw(:,:,:)
CALL tridia_solver( zwdia, zwsup, zwinf, ztw_in, ztw , 0 )
!
DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Interior value ( multiplied by wmask)
zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
......@@ -724,6 +730,8 @@ CONTAINS
REAL(wp) :: zwy_jm1, zfm_ui, zfm_ui_m1, zfm_vj, zfm_vj_m1, zfm_wk, zC2t_v, zC4t_v ! - -
REAL(wp) :: ztu, ztv, ztu_im1, ztu_ip1, ztv_jm1, ztv_jp1
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwx_3d, zwy_3d, zwz, ztw, zltu_3d, zltv_3d
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi_in, ztw_in ! Read-only copies to avoid INTENT conflicts
! in calls to tridia_solver
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz, zptry
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zwinf, zwdia, zwsup
LOGICAL :: ll_zAimp ! flag to apply adaptive implicit vertical advection
......@@ -739,6 +747,8 @@ CONTAINS
zwy_3d(:,:,:) = 0._wp
zwz(:,:,:) = 0._wp
zwi(:,:,:) = 0._wp
zwt(:,:,:) = 0._wp
!
l_trd = .FALSE. ! set local switches
l_hst = .FALSE.
......@@ -820,7 +830,10 @@ CONTAINS
END DO
IF ( ll_zAimp ) THEN
CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 )
! We need to use separate copies of zwi to avoid intent conflicts!
zwi_in(:,:,:) = zwi(:,:,:)
CALL tridia_solver( zwdia, zwsup, zwinf, zwi_in, zwi , 0 )
!
ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
DO_3D( 1, 1, 1, 1, 2, jpkm1 ) ! Interior value ( multiplied by wmask)
......@@ -934,7 +947,9 @@ CONTAINS
ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
END_3D
!
CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
! We need to use separate copies of ztw to avoid intent conflicts!
ztw_in(:,:,:) = ztw(:,:,:)
CALL tridia_solver( zwdia, zwsup, zwinf, ztw_in, ztw , 0 )
!
DO_3D( 1, 1, 1, 1, 2, jpkm1 ) ! Interior value ( multiplied by wmask)
zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
......
src/OCE/TRA/traqsr.F90
View file @ 672ff863
......@@ -48,6 +48,7 @@ MODULE traqsr
LOGICAL , PUBLIC :: ln_qsr_2bd !: 2 band light absorption flag
LOGICAL , PUBLIC :: ln_qsr_bio !: bio-model light absorption flag
INTEGER , PUBLIC :: nn_chldta !: use Chlorophyll data (=1) or not (=0)
REAL(wp), PUBLIC :: rn_chl_conc !: Chlorophyll concentration (for nn_chldta=0)
REAL(wp), PUBLIC :: rn_abs !: fraction absorbed in the very near surface (RGB & 2 bands)
REAL(wp), PUBLIC :: rn_si0 !: very near surface depth of extinction (RGB & 2 bands)
REAL(wp), PUBLIC :: rn_si1 !: deepest depth of extinction (water type I) (2 bands)
......@@ -221,7 +222,7 @@ CONTAINS
ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15
END_3D
ELSE !* constant chlorophyll
zchl = 0.05
zchl = rn_chl_conc
! NB. make sure constant value is such that:
zchl = MIN( 10. , MAX( 0.03, zchl ) )
! Convert chlorophyll value to attenuation coefficient look-up table index
......@@ -341,14 +342,14 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ios, irgb, ierror, ioptio ! local integer
REAL(wp) :: zz0, zc0 , zc1, zcoef ! local scalars
REAL(wp) :: zz1, zc2 , zc3, zchl ! - -
REAL(wp) :: zz0, zc0 , zc1, zcoef ! local scalars
REAL(wp) :: zz1, zc2 , zc3, zchl ! - -
!
CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files
TYPE(FLD_N) :: sn_chl ! informations about the chlorofyl field to be read
!!
NAMELIST/namtra_qsr/ sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, &
& nn_chldta, rn_abs, rn_si0, rn_si1
& nn_chldta, rn_chl_conc, rn_abs, rn_si0, rn_si1
!!----------------------------------------------------------------------
!
READ ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
......@@ -367,6 +368,7 @@ CONTAINS
WRITE(numout,*) ' 2 band light penetration ln_qsr_2bd = ', ln_qsr_2bd
WRITE(numout,*) ' bio-model light penetration ln_qsr_bio = ', ln_qsr_bio
WRITE(numout,*) ' RGB : Chl data (=1) or cst value (=0) nn_chldta = ', nn_chldta
WRITE(numout,*) ' Chlorophyll concentration (for nn_chldta=0) rn_chl_conc = ', rn_chl_conc
WRITE(numout,*) ' RGB & 2 bands: fraction of light (rn_si1) rn_abs = ', rn_abs
WRITE(numout,*) ' RGB & 2 bands: shortess depth of extinction rn_si0 = ', rn_si0
WRITE(numout,*) ' 2 bands: longest depth of extinction rn_si1 = ', rn_si1
......@@ -415,7 +417,7 @@ CONTAINS
& 'Solar penetration function of read chlorophyll', 'namtra_qsr' , no_print )
ENDIF
IF( nqsr == np_RGB ) THEN ! constant Chl
IF(lwp) WRITE(numout,*) ' ==>>> Constant Chlorophyll concentration = 0.05'
IF(lwp) WRITE(numout,*) ' ==>>> Constant Chlorophyll concentration = ', rn_chl_conc
ENDIF
!
CASE( np_2BD ) !== 2 bands light penetration ==!
......
src/OCE/ZDF/zdf_oce.F90
View file @ 672ff863
......@@ -39,6 +39,7 @@ MODULE zdf_oce
! ! gravity wave-induced vertical mixing
LOGICAL , PUBLIC :: ln_zdfswm !: surface wave-induced mixing flag
LOGICAL , PUBLIC :: ln_zdfiwm !: internal wave-induced mixing flag
LOGICAL , PUBLIC :: ln_zdftmx !: old tidal mixing scheme (Simmons et al 2004)
LOGICAL , PUBLIC :: ln_zdfmfc !: convection: eddy diffusivity Mass Flux Convection
! ! coefficients
REAL(wp), PUBLIC :: rn_avm0 !: vertical eddy viscosity (m2/s)
......
src/OCE/ZDF/zdfdrg.F90
View file @ 672ff863
......@@ -56,8 +56,9 @@ MODULE zdfdrg
LOGICAL :: ln_boost !: =T regional boost of Cd0 ; =F Cd0 horizontally uniform
REAL(wp) :: rn_boost !: local boost factor [ - ]
REAL(wp), PUBLIC :: r_Cdmin_top, r_Cdmax_top, r_z0_top, r_ke0_top ! set from namdrg_top namelist values
REAL(wp), PUBLIC :: r_Cdmin_bot, r_Cdmax_bot, r_z0_bot, r_ke0_bot ! - - namdrg_bot - -
LOGICAL :: ln_cdmin2d !: =T 2D varying Cdmin (ln_loglayer=T)
REAL(wp), PUBLIC :: r_Cdmax_top, r_z0_top, r_ke0_top ! set from namdrg_top namelist values
REAL(wp), PUBLIC :: r_Cdmax_bot, r_z0_bot, r_ke0_bot ! - - namdrg_bot - -
INTEGER :: ndrg ! choice of the type of drag coefficient
! ! associated indices:
......@@ -69,8 +70,9 @@ MODULE zdfdrg
LOGICAL , PUBLIC :: l_zdfdrg !: flag to update at each time step the top/bottom Cd
LOGICAL :: l_log_not_linssh !: flag to update at each time step the position ot the velocity point
!
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: rCd0_top, rCd0_bot !: precomputed top/bottom drag coeff. at t-point (>0)
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: rCdU_top, rCdU_bot !: top/bottom drag coeff. at t-point (<0) [m/s]
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: r_Cdmin_top, r_Cdmin_bot !: set from namelist values or file
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: rCd0_top, rCd0_bot !: precomputed top/bottom drag coeff. at t-point (>0)
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: rCdU_top, rCdU_bot !: top/bottom drag coeff. at t-point (<0) [m/s]
!! * Substitutions
# include "do_loop_substitute.h90"
......@@ -104,7 +106,7 @@ CONTAINS
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
! ! !! !== top or bottom variables ==!
INTEGER , DIMENSION(:,:), INTENT(in ) :: k_mk ! wet level (1st or last)
REAL(wp) , INTENT(in ) :: pCdmin ! min drag value
REAL(wp), DIMENSION(:,:), INTENT(in ) :: pCdmin ! min drag value
REAL(wp) , INTENT(in ) :: pCdmax ! max drag value
REAL(wp) , INTENT(in ) :: pz0 ! roughness
REAL(wp) , INTENT(in ) :: pke0 ! background tidal KE
......@@ -125,7 +127,7 @@ CONTAINS
!
!!JC: possible WAD implementation should modify line below if layers vanish
zcd = ( vkarmn / LOG( zzz / pz0 ) )**2
zcd = pCd0(ji,jj) * MIN( MAX( pCdmin , zcd ) , pCdmax ) ! here pCd0 = mask*boost
zcd = pCd0(ji,jj) * MIN( MAX( pCdmin(ji,jj) , zcd ) , pCdmax ) ! here pCd0 = mask*boost
pCdU(ji,jj) = - zcd * SQRT( 0.25 * ( zut*zut + zvt*zvt ) + pke0 )
END_2D
ELSE !== standard Cd ==!
......@@ -266,7 +268,7 @@ CONTAINS
!
! !== BOTTOM drag setting ==! (applied at seafloor)
!
ALLOCATE( rCd0_bot(jpi,jpj), rCdU_bot(jpi,jpj) )
ALLOCATE( rCd0_bot(jpi,jpj), rCdU_bot(jpi,jpj), r_Cdmin_bot(jpi,jpj) )
CALL drg_init( 'BOTTOM' , mbkt , & ! <== in
& r_Cdmin_bot, r_Cdmax_bot, r_z0_bot, r_ke0_bot, rCd0_bot, rCdU_bot ) ! ==> out
!
......@@ -277,7 +279,7 @@ CONTAINS
ENDIF
!
IF( ln_isfcav ) THEN
ALLOCATE( rCd0_top(jpi,jpj))
ALLOCATE( rCd0_top(jpi,jpj), r_Cdmin_top(jpi,jpj) )
CALL drg_init( 'TOP ' , mikt , & ! <== in
& r_Cdmin_top, r_Cdmax_top, r_z0_top, r_ke0_top, rCd0_top, rCdU_top ) ! ==> out
ENDIF
......@@ -295,21 +297,23 @@ CONTAINS
!!----------------------------------------------------------------------
CHARACTER(len=6) , INTENT(in ) :: cd_topbot ! top/ bot indicator
INTEGER , DIMENSION(:,:), INTENT(in ) :: k_mk ! 1st/last wet level
REAL(wp) , INTENT( out) :: pCdmin, pCdmax ! min and max drag coef. [-]
REAL(wp), DIMENSION(:,:), INTENT( out) :: pCdmin ! min drag coef. [-]
REAL(wp) , INTENT( out) :: pCdmax ! max drag coef. [-]
REAL(wp) , INTENT( out) :: pz0 ! roughness [m]
REAL(wp) , INTENT( out) :: pke0 ! background KE [m2/s2]
REAL(wp), DIMENSION(:,:), INTENT( out) :: pCd0 ! masked precomputed part of the non-linear drag coefficient
REAL(wp), DIMENSION(:,:), INTENT( out) :: pCdU ! minus linear drag*|U| at t-points [m/s]
!!
CHARACTER(len=40) :: cl_namdrg, cl_file, cl_varname, cl_namref, cl_namcfg ! local names
CHARACTER(len=40) :: cl_namdrg, cl_namref , cl_namcfg ! local names
CHARACTER(len=40) :: cl_file1 , cl_varname1, cl_file2 , cl_varname2 ! local names
INTEGER :: ji, jj ! dummy loop indexes
LOGICAL :: ll_top, ll_bot ! local logical
INTEGER :: ios, inum, imk ! local integers
REAL(wp):: zmsk, zzz, zcd ! local scalars
REAL(wp), DIMENSION(jpi,jpj) :: zmsk_boost ! 2D workspace
!!
NAMELIST/namdrg_top/ rn_Cd0, rn_Uc0, rn_Cdmax, rn_ke0, rn_z0, ln_boost, rn_boost
NAMELIST/namdrg_bot/ rn_Cd0, rn_Uc0, rn_Cdmax, rn_ke0, rn_z0, ln_boost, rn_boost
NAMELIST/namdrg_top/ rn_Cd0, rn_Uc0, rn_Cdmax, rn_ke0, rn_z0, ln_boost, rn_boost, ln_cdmin2d
NAMELIST/namdrg_bot/ rn_Cd0, rn_Uc0, rn_Cdmax, rn_ke0, rn_z0, ln_boost, rn_boost, ln_cdmin2d
!!----------------------------------------------------------------------
!
! !== set TOP / BOTTOM specificities ==!
......@@ -319,18 +323,22 @@ CONTAINS
SELECT CASE (cd_topbot)
CASE( 'TOP ' )
ll_top = .TRUE.
cl_namdrg = 'namdrg_top'
cl_namref = 'namdrg_top in reference namelist'
cl_namcfg = 'namdrg_top in configuration namelist'
cl_file = 'tfr_coef.nc'
cl_varname = 'tfr_coef'
cl_namdrg = 'namdrg_top'
cl_namref = 'namdrg_top in reference namelist'
cl_namcfg = 'namdrg_top in configuration namelist'
cl_file1 = 'tfr_coef.nc'
cl_varname1 = 'tfr_coef'
cl_file2 = 'tfr_cdmin_2d.nc'
cl_varname2 = 'tfr_cdmin'
CASE( 'BOTTOM' )
ll_bot = .TRUE.
cl_namdrg = 'namdrg_bot'
cl_namref = 'namdrg_bot in reference namelist'
cl_namcfg = 'namdrg_bot in configuration namelist'
cl_file = 'bfr_coef.nc'
cl_varname = 'bfr_coef'
cl_namdrg = 'namdrg_bot'
cl_namref = 'namdrg_bot in reference namelist'
cl_namcfg = 'namdrg_bot in configuration namelist'
cl_file1 = 'bfr_coef.nc'
cl_varname1 = 'bfr_coef'
cl_file2 = 'bfr_cdmin_2d.nc'
cl_varname2 = 'bfr_cdmin'
CASE DEFAULT
CALL ctl_stop( 'drg_init: bad value for cd_topbot ' )
END SELECT
......@@ -358,8 +366,17 @@ CONTAINS
WRITE(numout,*) ' associated boost factor rn_boost = ', rn_boost
ENDIF
!
! !== return some namelist parametres ==! (used in non_lin and loglayer cases)
pCdmin = rn_Cd0
! !== return some namelist parametres ==! (used in non_lin and loglayer cases)
!
IF( ln_cdmin2d ) THEN !* 2D varying pCdmin for the loglayer case
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' ==>>> use a 2D varying Cdmin read in ', TRIM(cl_file2), ' file'
CALL iom_open ( TRIM(cl_file2), inum )
CALL iom_get ( inum, jpdom_global, TRIM(cl_varname2), pCdmin, 1 )
CALL iom_close( inum)
ELSE !* no boost: boost factor = 1
pCdmin(:,:) = rn_Cd0
ENDIF
pCdmax = rn_Cdmax
pz0 = rn_z0
pke0 = rn_ke0
......@@ -368,11 +385,11 @@ CONTAINS
!
IF( ln_boost ) THEN !* regional boost: boost factor = 1 + regional boost
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' ==>>> use a regional boost read in ', TRIM(cl_file), ' file'
IF(lwp) WRITE(numout,*) ' ==>>> use a regional boost read in ', TRIM(cl_file1), ' file'
IF(lwp) WRITE(numout,*) ' using enhancement factor of ', rn_boost
! cl_varname is a coefficient in [0,1] giving where to apply the regional boost
CALL iom_open ( TRIM(cl_file), inum )
CALL iom_get ( inum, jpdom_global, TRIM(cl_varname), zmsk_boost, 1 )
CALL iom_open ( TRIM(cl_file1), inum )
CALL iom_get ( inum, jpdom_global, TRIM(cl_varname1), zmsk_boost, 1 )
CALL iom_close( inum)
zmsk_boost(:,:) = 1._wp + rn_boost * zmsk_boost(:,:)
!
......@@ -422,7 +439,13 @@ CONTAINS
IF(lwp) WRITE(numout,*) ' with a logarithmic Cd0 formulation Cd0 = ( vkarman log(z/z0) )^2 ,'
IF(lwp) WRITE(numout,*) ' a background velocity module of (rn_ke0)^1/2 = ', SQRT(pke0), 'm/s), '
IF(lwp) WRITE(numout,*) ' a logarithmic formulation: a roughness of ', pz0, ' meters, and '
IF(lwp) WRITE(numout,*) ' a proportionality factor bounded by min/max values of ', pCdmin, pCdmax
IF(lwp) WRITE(numout,*) ' a proportionality factor bounded by '
IF( ln_cdmin2d ) THEN
IF(lwp) WRITE(numout,*) ' a max value of ', pCdmax, 'and '
IF(lwp) WRITE(numout,*) ' a min 2D Cdmin scalar field read in ', TRIM(cl_file2), ' file'
ELSE
IF(lwp) WRITE(numout,*) ' min/max values of ', rn_Cd0, pCdmax
ENDIF
!
l_zdfdrg = .TRUE. !* Cd*|U| updated at each time-step (it depends on ocean velocity)
!
......
src/OCE/ZDF/zdfmxl.F90
View file @ 672ff863
......@@ -15,6 +15,7 @@ MODULE zdfmxl
USE dom_oce ! ocean space and time domain variables
USE trc_oce , ONLY: l_offline ! ocean space and time domain variables
USE zdf_oce ! ocean vertical physics
USE eosbn2 ! for zdf_mxl_zint
!
USE in_out_manager ! I/O manager
USE prtctl ! Print control
......@@ -26,15 +27,27 @@ MODULE zdfmxl
PRIVATE
PUBLIC zdf_mxl, zdf_mxl_turb, zdf_mxl_alloc ! called by zdfphy.F90
PUBLIC zdf_mxl_zint ! called by diahth.F90
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: nmln !: number of level in the mixed layer (used by LDF, ZDF, TRD, TOP)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmld !: mixing layer depth (turbocline) [m] (used by TOP)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlp !: mixed layer depth (rho=rho0+zdcrit) [m] (used by LDF)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlpt !: depth of the last T-point inside the mixed layer [m] (used by LDF)
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hmld_zint !: vertically-interpolated mixed layer depth [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: htc_mld ! Heat content of hmld_zint
LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ll_found ! Is T_b to be found by interpolation ?
LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ll_belowml ! Flag points below mixed layer when ll_found=F
REAL(wp), PUBLIC :: rho_c = 0.01_wp !: density criterion for mixed layer depth
REAL(wp), PUBLIC :: avt_c = 5.e-4_wp ! Kz criterion for the turbocline depth
TYPE, PUBLIC :: MXL_ZINT !: Structure for MLD defs
INTEGER :: mld_type ! mixed layer type
REAL(wp) :: zref ! depth of initial T_ref
REAL(wp) :: dT_crit ! Critical temp diff
REAL(wp) :: iso_frac ! Fraction of rn_dT_crit
END TYPE MXL_ZINT
!! * Substitutions
# include "do_loop_substitute.h90"
# include "single_precision_substitute.h90"
......@@ -52,7 +65,8 @@ CONTAINS
!!----------------------------------------------------------------------
zdf_mxl_alloc = 0 ! set to zero if no array to be allocated
IF( .NOT. ALLOCATED( nmln ) ) THEN
ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc )
ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj), &
& htc_mld(jpi,jpj), ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc )
!
CALL mpp_sum ( 'zdfmxl', zdf_mxl_alloc )
IF( zdf_mxl_alloc /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl_alloc: failed to allocate arrays.' )
......@@ -86,6 +100,8 @@ CONTAINS
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth'
IF(lwp) WRITE(numout,*) '~~~~~~~ '
! ! allocate zdfmxl arrays
IF( zdf_mxl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl : unable to allocate arrays' )
ENDIF
ENDIF
!
......@@ -120,6 +136,352 @@ CONTAINS
!
END SUBROUTINE zdf_mxl
SUBROUTINE zdf_mxl_zint_mld( sf , Kmm)
!!----------------------------------------------------------------------------------
!! *** ROUTINE zdf_mxl_zint_mld ***
!
! Calculate vertically-interpolated mixed layer depth diagnostic.
!
! This routine can calculate the mixed layer depth diagnostic suggested by
! Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate
! vertically-interpolated mixed-layer depth diagnostics with other parameter
! settings set in the namzdf_mldzint namelist.
!
! If mld_type=1 the mixed layer depth is calculated as the depth at which the
! density has increased by an amount equivalent to a temperature difference of
! 0.8C at the surface.
!
! For other values of mld_type the mixed layer is calculated as the depth at
! which the temperature differs by 0.8C from the surface temperature.
!
! David Acreman, Daley Calvert
!
!!-----------------------------------------------------------------------------------
TYPE(MXL_ZINT), INTENT(in) :: sf
INTEGER, INTENT(in) :: Kmm ! ocean time level index
! Diagnostic criteria
INTEGER :: nn_mld_type ! mixed layer type
REAL(wp) :: rn_zref ! depth of initial T_ref
REAL(wp) :: rn_dT_crit ! Critical temp diff
REAL(wp) :: rn_iso_frac ! Fraction of rn_dT_crit used
! Local variables
REAL(wp), PARAMETER :: zepsilon = 1.e-30 ! local small value
INTEGER, DIMENSION(jpi,jpj) :: ikmt ! number of active tracer levels
INTEGER, DIMENSION(jpi,jpj) :: ik_ref ! index of reference level
INTEGER, DIMENSION(jpi,jpj) :: ik_iso ! index of last uniform temp level
REAL, DIMENSION(jpi,jpj,jpk) :: zT ! Temperature or density
REAL, DIMENSION(jpi,jpj) :: ppzdep ! depth for use in calculating d(rho)
REAL, DIMENSION(jpi,jpj) :: zT_ref ! reference temperature
REAL :: zT_b ! base temperature
REAL, DIMENSION(jpi,jpj,jpk) :: zdTdz ! gradient of zT
REAL, DIMENSION(jpi,jpj,jpk) :: zmoddT ! Absolute temperature difference
REAL :: zdz ! depth difference
REAL :: zdT ! temperature difference
REAL, DIMENSION(jpi,jpj) :: zdelta_T ! difference critereon
REAL, DIMENSION(jpi,jpj) :: zRHO1, zRHO2 ! Densities
INTEGER :: ji, jj, jk ! loop counter
!!-------------------------------------------------------------------------------------
!
! Unpack structure
nn_mld_type = sf%mld_type
rn_zref = sf%zref
rn_dT_crit = sf%dT_crit
rn_iso_frac = sf%iso_frac
! Set the mixed layer depth criterion at each grid point
IF( nn_mld_type == 0 ) THEN
zdelta_T(:,:) = rn_dT_crit
zT(:,:,:) = rhop(:,:,:)
ELSE IF( nn_mld_type == 1 ) THEN
ppzdep(:,:)=0.0
call eos ( ts(:,:,1,:,Kmm), ppzdep(:,:), zRHO1(:,:) )
! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST
! [assumes number of tracers less than number of vertical levels]
zT(:,:,1:jpts)=ts(:,:,1,1:jpts,Kmm)
zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit
CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) )
zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rho0
! RHO from eos (2d version) doesn't calculate north or east halo:
CALL lbc_lnk( 'zdfmxl', zdelta_T, 'T', 1. )
zT(:,:,:) = rhop(:,:,:)
ELSE
zdelta_T(:,:) = rn_dT_crit
zT(:,:,:) = ts(:,:,:,jp_tem,Kmm)
END IF
! Calculate the gradient of zT and absolute difference for use later
DO jk = 1 ,jpk-2
zdTdz(:,:,jk) = ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w(:,:,jk+1,Kmm)
zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )
END DO
! Find density/temperature at the reference level (Kara et al use 10m).
! ik_ref is the index of the box centre immediately above or at the reference level
! Find rn_zref in the array of model level depths and find the ref
! density/temperature by linear interpolation.
DO jk = jpkm1, 2, -1
WHERE ( gdept(:,:,jk,Kmm) > rn_zref )
ik_ref(:,:) = jk - 1
zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - gdept(:,:,jk-1,Kmm) )
END WHERE
END DO
! If the first grid box centre is below the reference level then use the
! top model level to get zT_ref
WHERE ( gdept(:,:,1,Kmm) > rn_zref )
zT_ref = zT(:,:,1)
ik_ref = 1
END WHERE
! The number of active tracer levels is 1 less than the number of active w levels
ikmt(:,:) = mbkt(:,:) - 1
! Initialize / reset
ll_found(:,:) = .false.
IF ( rn_iso_frac - zepsilon > 0. ) THEN
! Search for a uniform density/temperature region where adjacent levels
! differ by less than rn_iso_frac * deltaT.
! ik_iso is the index of the last level in the uniform layer
! ll_found indicates whether the mixed layer depth can be found by interpolation
ik_iso(:,:) = ik_ref(:,:)
DO jj = 1, jpj ! Changed from nlcj
DO ji = 1, jpi ! Changed from nlci
!CDIR NOVECTOR
DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1
IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN
ik_iso(ji,jj) = jk
ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )
EXIT
END IF
END DO
END DO
END DO
! Use linear interpolation to find depth of mixed layer base where possible
hmld_zint(:,:) = rn_zref
DO jj = 1, jpj
DO ji = 1, jpi
IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN
zdz = abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )
hmld_zint(ji,jj) = gdept(ji,jj,ik_iso(ji,jj),Kmm) + zdz
END IF
END DO
END DO
END IF
! If ll_found = .false. then calculate MLD using difference of zdelta_T
! from the reference density/temperature
! Prevent this section from working on land points
WHERE ( tmask(:,:,1) /= 1.0 )
ll_found = .true.
END WHERE
DO jk=1, jpk
ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)
END DO
! Set default value where interpolation cannot be used (ll_found=false)
DO jj = 1, jpj
DO ji = 1, jpi
IF ( .not. ll_found(ji,jj) ) hmld_zint(ji,jj) = gdept(ji,jj,ikmt(ji,jj),Kmm)
END DO
END DO
DO jj = 1, jpj
DO ji = 1, jpi
!CDIR NOVECTOR
DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj)
IF ( ll_found(ji,jj) ) EXIT
IF ( ll_belowml(ji,jj,jk) ) THEN
zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) )
zdT = zT_b - zT(ji,jj,jk-1)
zdz = zdT / zdTdz(ji,jj,jk-1)
hmld_zint(ji,jj) = gdept(ji,jj,jk-1,Kmm) + zdz
EXIT
END IF
END DO
END DO
END DO
hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1)
!
END SUBROUTINE zdf_mxl_zint_mld
SUBROUTINE zdf_mxl_zint_htc( kt , Kmm)
!!----------------------------------------------------------------------
!! *** ROUTINE zdf_mxl_zint_htc ***
!!
!! ** Purpose :
!!
!! ** Method :
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time-step index
INTEGER, INTENT(in) :: Kmm ! ocean time level index
INTEGER :: ji, jj, jk
INTEGER :: ikmax
REAL(wp) :: zc, zcoef
!
INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ilevel
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zthick_0, zthick
!!----------------------------------------------------------------------
IF( .NOT. ALLOCATED(ilevel) ) THEN
ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), &
& zthick(jpi,jpj), STAT=ji )
IF( lk_mpp ) CALL mpp_sum( 'zdfmxl', ji )
IF( ji /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' )
ENDIF
! Find last whole model T level above the MLD
ilevel(:,:) = 0
zthick_0(:,:) = 0._wp
DO jk = 1, jpkm1
DO jj = 1, jpj
DO ji = 1, jpi
zthick_0(ji,jj) = zthick_0(ji,jj) + e3t(ji,jj,jk,Kmm)
IF( zthick_0(ji,jj) < hmld_zint(ji,jj) ) ilevel(ji,jj) = jk
END DO
END DO
WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2)
WRITE(numout,*) 'gdepw(jk+1 =',jk+1,') =',gdepw(2,2,jk+1,Kmm)
END DO
! Surface boundary condition
IF( ln_linssh ) THEN ; zthick(:,:) = ssh(:,:,Kmm) ; htc_mld(:,:) = ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) * tmask(:,:,1)
ELSE ; zthick(:,:) = 0._wp ; htc_mld(:,:) = 0._wp
ENDIF
! Deepest whole T level above the MLD
ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 )
! Integration down to last whole model T level
DO jk = 1, ikmax
DO jj = 1, jpj
DO ji = 1, jpi
zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1 ) ) ! 0 below ilevel
zthick(ji,jj) = zthick(ji,jj) + zc
htc_mld(ji,jj) = htc_mld(ji,jj) + zc * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
END DO
END DO
END DO
! Subsequent partial T level
zthick(:,:) = hmld_zint(:,:) - zthick(:,:) ! remaining thickness to reach MLD
DO jj = 1, jpj
DO ji = 1, jpi
htc_mld(ji,jj) = htc_mld(ji,jj) + ts(ji,jj,ilevel(ji,jj)+1,jp_tem,Kmm) &
& * MIN( e3t(ji,jj,ilevel(ji,jj)+1,Kmm), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1)
END DO
END DO
WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2)
! Convert to heat content
zcoef = rho0 * rcp
htc_mld(:,:) = zcoef * htc_mld(:,:)
END SUBROUTINE zdf_mxl_zint_htc
SUBROUTINE zdf_mxl_zint( kt , Kmm)
!!----------------------------------------------------------------------
!! *** ROUTINE zdf_mxl_zint ***
!!
!! ** Purpose :
!!
!! ** Method :
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time-step index
INTEGER, INTENT(in) :: Kmm ! ocean time level index
INTEGER :: ios
INTEGER :: jn
INTEGER :: nn_mld_diag = 0 ! number of diagnostics
CHARACTER(len=1) :: cmld
LOGICAL, SAVE, DIMENSION(5) :: l_mld, l_htc
TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5
TYPE(MXL_ZINT), SAVE, DIMENSION(5) :: mld_diags
NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5
!!----------------------------------------------------------------------
IF( kt == nit000 ) THEN
READ ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist' )
READ ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 )
902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist' )
IF(lwm) WRITE ( numond, namzdf_mldzint )
IF( nn_mld_diag > 5 ) CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' )
l_mld(:) = .FALSE. ; l_htc(:) = .FALSE.
mld_diags(1) = sn_mld1
mld_diags(2) = sn_mld2
mld_diags(3) = sn_mld3
mld_diags(4) = sn_mld4
mld_diags(5) = sn_mld5
IF( nn_mld_diag > 0 ) THEN
IF( lwp ) THEN
WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================'
WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)'
ENDIF
DO jn = 1, nn_mld_diag
! Whether the diagnostic is requested
WRITE(cmld,'(I1)') jn
IF( iom_use( "mldzint_"//cmld ) ) l_mld(jn) = .TRUE.
IF( iom_use( "mldhtc_"//cmld ) ) l_htc(jn) = .TRUE.
IF( lwp ) THEN
WRITE(numout,*) 'MLD criterion',jn,':'
WRITE(numout,*) ' nn_mld_type =', mld_diags(jn)%mld_type
WRITE(numout,*) ' rn_zref =' , mld_diags(jn)%zref
WRITE(numout,*) ' rn_dT_crit =' , mld_diags(jn)%dT_crit
WRITE(numout,*) ' rn_iso_frac =', mld_diags(jn)%iso_frac
ENDIF
END DO
WRITE(numout,*) '===================================================================='
ENDIF
ENDIF
IF( nn_mld_diag > 0 ) THEN
DO jn = 1, nn_mld_diag
WRITE(cmld,'(I1)') jn
IF( l_mld(jn) .OR. l_htc(jn) ) THEN
CALL zdf_mxl_zint_mld( mld_diags(jn), Kmm)
IF( l_mld(jn) ) CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) )
IF( l_htc(jn) ) THEN
CALL zdf_mxl_zint_htc( kt, Kmm )
CALL iom_put( "mldhtc_"//cmld, htc_mld(:,:) )
ENDIF
ENDIF
END DO
ENDIF
END SUBROUTINE zdf_mxl_zint
SUBROUTINE zdf_mxl_turb( kt, Kmm )
!!----------------------------------------------------------------------
......