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
  • nemo/nemo
  • sparonuz/nemo
  • hatfield/nemo
  • extdevs/nemo
4 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 657 additions and 710 deletions
src/ICE/iceupdate.F90
View file @ 78481193
......@@ -314,7 +314,7 @@ CONTAINS
!!
!! ** Action : * at each ice time step (every nn_fsbc time step):
!! - compute the modulus of ice-ocean relative velocity
!! (*rho*Cd) at T-point (C-grid) or I-point (B-grid)
!! (*rho*Cd) at T-point (C-grid)
!! tmod_io = rhoco * | U_ice-U_oce |
!! - update the modulus of stress at ocean surface
!! taum = (1-a) * taum + a * tmod_io * | U_ice-U_oce |
......@@ -325,19 +325,19 @@ CONTAINS
!!
!! NB: - ice-ocean rotation angle no more allowed
!! - here we make an approximation: taum is only computed every ice time step
!! This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
!! to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
!! This avoids mutiple average to pass from U,V grids to T grids
!! taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
!!
!! ** Outputs : - utau, vtau : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
!! ** Outputs : - utau, vtau : surface ocean i- and j-stress (T-pts) updated with ice-ocean fluxes
!! - taum : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
!!---------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! ocean time-step index
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pu_oce, pv_oce ! surface ocean currents
!
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zat_u, zutau_ice, zu_t, zmodt ! local scalar
REAL(wp) :: zat_v, zvtau_ice, zv_t, zrhoco ! - -
REAL(wp) :: zflagi ! - -
REAL(wp) :: zutau_ice, zu_t, zmodt ! local scalar
REAL(wp) :: zvtau_ice, zv_t, zrhoco ! - -
REAL(wp) :: zflagi ! - -
!!---------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('iceupdate')
......@@ -352,8 +352,8 @@ CONTAINS
IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN !== Ice time-step only ==! (i.e. surface module time-step)
DO_2D( 0, 0, 0, 0 ) !* update the modulus of stress at ocean surface (T-point)
! ! 2*(U_ice-U_oce) at T-point
zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)
zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1)
zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj) ! u_oce = ssu_m
zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1) ! v_oce = ssv_m
! ! |U_ice-U_oce|^2
zmodt = 0.25_wp * ( zu_t * zu_t + zv_t * zv_t )
! ! update the ocean stress modulus
......@@ -377,19 +377,14 @@ CONTAINS
ENDIF
!
DO_2D( 0, 0, 0, 0 ) !* update the stress WITHOUT an ice-ocean rotation angle
! ice area at u and v-points
zat_u = ( at_i(ji,jj) * tmask(ji,jj,1) + at_i (ji+1,jj ) * tmask(ji+1,jj ,1) ) &
& / MAX( 1.0_wp , tmask(ji,jj,1) + tmask(ji+1,jj ,1) )
zat_v = ( at_i(ji,jj) * tmask(ji,jj,1) + at_i (ji ,jj+1 ) * tmask(ji ,jj+1,1) ) &
& / MAX( 1.0_wp , tmask(ji,jj,1) + tmask(ji ,jj+1,1) )
! ! linearized quadratic drag formulation
zutau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
zvtau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
zutau_ice = 0.5_wp * tmod_io(ji,jj) * ( u_ice(ji,jj) + u_ice(ji-1,jj) - pu_oce(ji,jj) - pu_oce(ji-1,jj) )
zvtau_ice = 0.5_wp * tmod_io(ji,jj) * ( v_ice(ji,jj) + v_ice(ji,jj-1) - pv_oce(ji,jj) - pv_oce(ji,jj-1) )
! ! stresses at the ocean surface
utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
utau(ji,jj) = ( 1._wp - at_i(ji,jj) ) * utau_oce(ji,jj) + at_i(ji,jj) * zutau_ice
vtau(ji,jj) = ( 1._wp - at_i(ji,jj) ) * vtau_oce(ji,jj) + at_i(ji,jj) * zvtau_ice
END_2D
CALL lbc_lnk( 'iceupdate', utau, 'U', -1.0_wp, vtau, 'V', -1.0_wp ) ! lateral boundary condition
CALL lbc_lnk( 'iceupdate', utau, 'T', -1.0_wp, vtau, 'T', -1.0_wp ) ! lateral boundary condition
!
IF( ln_timing ) CALL timing_stop('iceupdate')
!
......
src/OCE/CRS/crsfld.F90
View file @ 78481193
......@@ -211,8 +211,8 @@ CONTAINS
! sbc fields
CALL crs_dom_ope( ssh(:,:,Kmm) , 'VOL', 'T', tmask, sshn_crs , p_e12=e1e2t, p_e3=ze3t , psgn=1.0_wp )
CALL crs_dom_ope( utau , 'SUM', 'U', umask, utau_crs , p_e12=e2u , p_surf_crs=e2u_crs , psgn=1.0_wp )
CALL crs_dom_ope( vtau , 'SUM', 'V', vmask, vtau_crs , p_e12=e1v , p_surf_crs=e1v_crs , psgn=1.0_wp )
CALL crs_dom_ope( utau , 'SUM', 'T', tmask, utau_crs , p_e12=e2t , p_surf_crs=e2t_crs , psgn=1.0_wp ) !clem tau: check psgn ??
CALL crs_dom_ope( vtau , 'SUM', 'T', tmask, vtau_crs , p_e12=e1t , p_surf_crs=e1t_crs , psgn=1.0_wp ) !
CALL crs_dom_ope( wndm , 'SUM', 'T', tmask, wndm_crs , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0_wp )
CALL crs_dom_ope( rnf , 'MAX', 'T', tmask, rnf_crs , psgn=1.0_wp )
CALL crs_dom_ope( qsr , 'SUM', 'T', tmask, qsr_crs , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0_wp )
......
src/OCE/DIA/diadetide.F90
View file @ 78481193
......@@ -5,11 +5,11 @@ MODULE diadetide
!!======================================================================
!! History : ! 2019 (S. Mueller)
!!----------------------------------------------------------------------
USE par_oce , ONLY : wp, jpi, jpj
USE in_out_manager , ONLY : lwp, numout
USE iom , ONLY : iom_put
USE dom_oce , ONLY : rn_Dt, nsec_day
USE phycst , ONLY : rpi
USE par_oce
USE in_out_manager
USE iom
USE dom_oce
USE phycst
USE tide_mod
#if defined key_xios
USE xios
......@@ -24,6 +24,8 @@ MODULE diadetide
PUBLIC :: dia_detide_init, dia_detide
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2019)
!! $Id$
......@@ -90,9 +92,9 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt
REAL(wp), DIMENSION(jpi,jpj) :: zwght_2D
REAL(wp), DIMENSION(A2D(0)) :: zwght_2D
REAL(wp) :: zwght, ztmp
INTEGER :: jn
INTEGER :: ji, jj, jn
! Compute detiding weight at the current time-step; the daily total weight
! is one, and the daily summation of a diagnosed field multiplied by this
......@@ -104,7 +106,10 @@ CONTAINS
zwght = zwght + 1.0_wp / REAL( ndiadetide, KIND=wp )
END IF
END DO
zwght_2D(:,:) = zwght
DO_2D( 0, 0, 0, 0 )
zwght_2D(ji,jj) = zwght
END_2D
CALL iom_put( "diadetide_weight", zwght_2D)
END SUBROUTINE dia_detide
......
src/OCE/DIA/diahsb.F90
View file @ 78481193
......@@ -50,6 +50,7 @@ MODULE diahsb
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmask_ini
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
......@@ -82,30 +83,61 @@ CONTAINS
REAL(wp) :: z_frc_trd_v ! - -
REAL(wp) :: z_wn_trd_t , z_wn_trd_s ! - -
REAL(wp) :: z_ssh_hc , z_ssh_sc ! - -
REAL(wp), DIMENSION(jpi,jpj,13) :: ztmp
REAL(wp), DIMENSION(jpi,jpj,jpkm1,4) :: ztmpk
REAL(wp), DIMENSION(A2D(0),13) :: ztmp
REAL(wp), DIMENSION(A2D(0),jpkm1,4) :: ztmpk
REAL(wp), DIMENSION(17) :: zbg
!!---------------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('dia_hsb')
!
ztmp (:,:,:) = 0._wp ! should be better coded
ztmpk(:,:,:,:) = 0._wp ! should be better coded
!
ts(:,:,:,1,Kmm) = ts(:,:,:,1,Kmm) * tmask(:,:,:) ; ts(:,:,:,1,Kbb) = ts(:,:,:,1,Kbb) * tmask(:,:,:) ;
ts(:,:,:,2,Kmm) = ts(:,:,:,2,Kmm) * tmask(:,:,:) ; ts(:,:,:,2,Kbb) = ts(:,:,:,2,Kbb) * tmask(:,:,:) ;
DO_2D( 0, 0, 0, 0 )
ztmp (ji,jj,:) = 0._wp ! should be better coded
ztmpk(ji,jj,:,:) = 0._wp ! should be better coded
!
ts(ji,jj,:,1,Kmm) = ts(ji,jj,:,1,Kmm) * tmask(ji,jj,:)
ts(ji,jj,:,1,Kbb) = ts(ji,jj,:,1,Kbb) * tmask(ji,jj,:)
!
ts(ji,jj,:,2,Kmm) = ts(ji,jj,:,2,Kmm) * tmask(ji,jj,:)
ts(ji,jj,:,2,Kbb) = ts(ji,jj,:,2,Kbb) * tmask(ji,jj,:)
END_2D
!
! ------------------------- !
! 1 - Trends due to forcing !
! ------------------------- !
! prepare trends
ztmp(:,:,1) = - r1_rho0 * ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) ) * surf(:,:) ! volume
ztmp(:,:,2) = sbc_tsc(:,:,jp_tem) * surf(:,:) ! heat
ztmp(:,:,3) = sbc_tsc(:,:,jp_sal) * surf(:,:) ! salt
IF( ln_rnf ) ztmp(:,:,4) = rnf_tsc(:,:,jp_tem) * surf(:,:) ! runoff temp
IF( ln_rnf_sal ) ztmp(:,:,5) = rnf_tsc(:,:,jp_sal) * surf(:,:) ! runoff salt
IF( ln_isf ) ztmp(:,:,6) = ( risf_cav_tsc(:,:,jp_tem) + risf_par_tsc(:,:,jp_tem) ) * surf(:,:) ! isf temp
IF( ln_traqsr ) ztmp(:,:,7) = r1_rho0_rcp * qsr(:,:) * surf(:,:) ! penetrative solar radiation
IF( ln_trabbc ) ztmp(:,:,8) = qgh_trd0(:,:) * surf(:,:) ! geothermal heat
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,1) = - r1_rho0 * ( emp(ji,jj) & ! volume
& - rnf(ji,jj) &
& - fwfisf_cav(ji,jj) &
& - fwfisf_par(ji,jj) ) * surf(ji,jj)
ztmp(ji,jj,2) = sbc_tsc(ji,jj,jp_tem) * surf(ji,jj) ! heat
ztmp(ji,jj,3) = sbc_tsc(ji,jj,jp_sal) * surf(ji,jj) ! salt
END_2D
IF( ln_rnf ) THEN
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,4) = rnf_tsc(ji,jj,jp_tem) * surf(ji,jj) ! runoff temp
END_2D
END IF
IF( ln_rnf_sal ) THEN
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,5) = rnf_tsc(ji,jj,jp_sal) * surf(ji,jj) ! runoff salt
END_2D
END IF
IF( ln_isf ) THEN
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,6) = ( risf_cav_tsc(ji,jj,jp_tem) &
& + risf_par_tsc(ji,jj,jp_tem) ) * surf(ji,jj) ! isf temp
END_2D
END IF
IF( ln_traqsr ) THEN
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,7) = r1_rho0_rcp * qsr(ji,jj) * surf(ji,jj) ! penetrative solar radiation
END_2D
END IF
IF( ln_trabbc ) THEN
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,8) = qgh_trd0(ji,jj) * surf(ji,jj) ! geothermal heat
END_2D
END IF
!
IF( ln_linssh ) THEN ! Advection flux through fixed surface (z=0)
IF( ln_isfcav ) THEN
......@@ -116,8 +148,10 @@ CONTAINS
END DO
END DO
ELSE
ztmp(:,:,9 ) = - surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_tem,Kbb)
ztmp(:,:,10) = - surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_sal,Kbb)
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,9 ) = - surf(ji,jj) * ww(ji,jj,1) * ts(ji,jj,1,jp_tem,Kbb)
ztmp(ji,jj,10) = - surf(ji,jj) * ww(ji,jj,1) * ts(ji,jj,1,jp_sal,Kbb)
END_2D
END IF
ENDIF
......@@ -152,7 +186,9 @@ CONTAINS
! glob_sum is needed because you keep only the interior domain to compute the sum (iscpl)
!
! ! volume variation (calculated with ssh)
ztmp(:,:,11) = surf(:,:)*ssh(:,:,Kmm) - surf_ini(:,:)*ssh_ini(:,:)
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,11) = surf(ji,jj)*ssh(ji,jj,Kmm) - surf_ini(ji,jj)*ssh_ini(ji,jj)
END_2D
! ! heat & salt content variation (associated with ssh)
IF( ln_linssh ) THEN ! linear free surface case
......@@ -164,8 +200,10 @@ CONTAINS
END DO
END DO
ELSE ! no under ice-shelf seas
ztmp(:,:,12) = surf(:,:) * ( ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) - ssh_hc_loc_ini(:,:) )
ztmp(:,:,13) = surf(:,:) * ( ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) - ssh_sc_loc_ini(:,:) )
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,12) = surf(ji,jj) * ( ts(ji,jj,1,jp_tem,Kmm) * ssh(ji,jj,Kmm) - ssh_hc_loc_ini(ji,jj) )
ztmp(ji,jj,13) = surf(ji,jj) * ( ts(ji,jj,1,jp_sal,Kmm) * ssh(ji,jj,Kmm) - ssh_sc_loc_ini(ji,jj) )
END_2D
END IF
ENDIF
......@@ -185,19 +223,27 @@ CONTAINS
! glob_sum is needed because you keep only the interior domain to compute the sum (iscpl)
!
DO jk = 1, jpkm1 ! volume
ztmpk(:,:,jk,1) = surf (:,:) * e3t(:,:,jk,Kmm)*tmask(:,:,jk) &
& - surf_ini(:,:) * e3t_ini(:,:,jk )*tmask_ini(:,:,jk)
DO_2D( 0, 0, 0, 0 )
ztmpk(ji,jj,jk,1) = surf (ji,jj) * e3t(ji,jj,jk,Kmm)*tmask(ji,jj,jk) &
& - surf_ini(ji,jj) * e3t_ini(ji,jj,jk )*tmask_ini(ji,jj,jk)
END_2D
END DO
DO jk = 1, jpkm1 ! heat
ztmpk(:,:,jk,2) = ( surf (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_tem,Kmm) &
& - surf_ini(:,:) * hc_loc_ini(:,:,jk) )
DO_2D( 0, 0, 0, 0 )
ztmpk(ji,jj,jk,2) = ( surf (ji,jj) * e3t(ji,jj,jk,Kmm)*ts(ji,jj,jk,jp_tem,Kmm) &
& - surf_ini(ji,jj) * hc_loc_ini(ji,jj,jk) )
END_2D
END DO
DO jk = 1, jpkm1 ! salt
ztmpk(:,:,jk,3) = ( surf (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_sal,Kmm) &
& - surf_ini(:,:) * sc_loc_ini(:,:,jk) )
DO_2D( 0, 0, 0, 0 )
ztmpk(ji,jj,jk,3) = ( surf (ji,jj) * e3t(ji,jj,jk,Kmm)*ts(ji,jj,jk,jp_sal,Kmm) &
& - surf_ini(ji,jj) * sc_loc_ini(ji,jj,jk) )
END_2D
END DO
DO jk = 1, jpkm1 ! total ocean volume
ztmpk(:,:,jk,4) = surf(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
DO_2D( 0, 0, 0, 0 )
ztmpk(ji,jj,jk,4) = surf(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
END_2D
END DO
! global sum
......@@ -315,14 +361,18 @@ CONTAINS
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' dia_hsb_rst : initialise hsb at initial state '
IF(lwp) WRITE(numout,*)
surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:) ! initial ocean surface
ssh_ini(:,:) = ssh(:,:,Kmm) ! initial ssh
DO_2D( 0, 0, 0, 0 )
surf_ini(ji,jj) = e1e2t(ji,jj) * tmask_i(ji,jj) ! initial ocean surface
ssh_ini(ji,jj) = ssh(ji,jj,Kmm) ! initial ssh
END_2D
DO jk = 1, jpk
! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance).
e3t_ini (:,:,jk) = e3t(:,:,jk,Kmm) * tmask(:,:,jk) ! initial vertical scale factors
tmask_ini (:,:,jk) = tmask(:,:,jk) ! initial mask
hc_loc_ini(:,:,jk) = ts(:,:,jk,jp_tem,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) ! initial heat content
sc_loc_ini(:,:,jk) = ts(:,:,jk,jp_sal,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) ! initial salt content
DO_2D( 0, 0, 0, 0 )
e3t_ini (ji,jj,jk) = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) ! initial vertical scale factors
tmask_ini (ji,jj,jk) = tmask(ji,jj,jk) ! initial mask
hc_loc_ini(ji,jj,jk) = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) ! initial heat content
sc_loc_ini(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) ! initial salt content
END_2D
END DO
frc_v = 0._wp ! volume trend due to forcing
frc_t = 0._wp ! heat content - - - -
......@@ -334,13 +384,15 @@ CONTAINS
ssh_hc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) ! initial heat content in ssh
ssh_sc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) ! initial salt content in ssh
END DO
END DO
ELSE
ssh_hc_loc_ini(:,:) = ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) ! initial heat content in ssh
ssh_sc_loc_ini(:,:) = ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) ! initial salt content in ssh
END DO
ELSE
DO_2D( 0, 0, 0, 0 )
ssh_hc_loc_ini(ji,jj) = ts(ji,jj,1,jp_tem,Kmm) * ssh(ji,jj,Kmm) ! initial heat content in ssh
ssh_sc_loc_ini(ji,jj) = ts(ji,jj,1,jp_sal,Kmm) * ssh(ji,jj,Kmm) ! initial salt content in ssh
END_2D
END IF
frc_wn_t = 0._wp ! initial heat content misfit due to free surface
frc_wn_s = 0._wp ! initial salt content misfit due to free surface
frc_wn_t = 0._wp ! initial heat content misfit due to free surface
frc_wn_s = 0._wp ! initial salt content misfit due to free surface
ENDIF
ENDIF
!
......@@ -388,6 +440,7 @@ CONTAINS
INTEGER, INTENT(in) :: Kmm ! time level index
!
INTEGER :: ierror, ios ! local integer
INTEGER :: ji, jj ! loop index
!!
NAMELIST/namhsb/ ln_diahsb
!!----------------------------------------------------------------------
......@@ -427,7 +480,10 @@ CONTAINS
! ----------------------------------------------- !
! 2 - Time independant variables and file opening !
! ----------------------------------------------- !
surf(:,:) = e1e2t(:,:) * tmask_i(:,:) ! masked surface grid cell area
DO_2D( 0, 0, 0, 0 )
surf(ji,jj) = e1e2t(ji,jj) * tmask_i(ji,jj) ! masked surface grid cell area
END_2D
surf_tot = glob_sum( 'diahsb', surf(:,:) ) ! total ocean surface area
IF( ln_bdy ) CALL ctl_warn( 'dia_hsb_init: heat/salt budget does not consider open boundary fluxes' )
......
src/OCE/DIA/diahth.F90
View file @ 78481193
......@@ -86,22 +86,22 @@ CONTAINS
INTEGER, INTENT( in ) :: kt ! ocean time-step index
INTEGER, INTENT( in ) :: Kmm ! ocean time level index
!!
INTEGER :: ji, jj, jk ! dummy loop arguments
REAL(wp) :: zrho3 = 0.03_wp ! density criterion for mixed layer depth
REAL(wp) :: zrho1 = 0.01_wp ! density criterion for mixed layer depth
REAL(wp) :: ztem2 = 0.2_wp ! temperature criterion for mixed layer depth
REAL(wp) :: zztmp, zzdep ! temporary scalars inside do loop
REAL(wp) :: zu, zv, zw, zut, zvt ! temporary workspace
REAL(wp), DIMENSION(jpi,jpj) :: zabs2 ! MLD: abs( tn - tn(10m) ) = ztem2
REAL(wp), DIMENSION(jpi,jpj) :: ztm2 ! Top of thermocline: tn = tn(10m) - ztem2
REAL(wp), DIMENSION(jpi,jpj) :: zrho10_3 ! MLD: rho = rho10m + zrho3
REAL(wp), DIMENSION(jpi,jpj) :: zpycn ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
REAL(wp), DIMENSION(jpi,jpj) :: ztinv ! max of temperature inversion
REAL(wp), DIMENSION(jpi,jpj) :: zdepinv ! depth of temperature inversion
REAL(wp), DIMENSION(jpi,jpj) :: zrho0_3 ! MLD rho = rho(surf) = 0.03
REAL(wp), DIMENSION(jpi,jpj) :: zrho0_1 ! MLD rho = rho(surf) = 0.01
REAL(wp), DIMENSION(jpi,jpj) :: zmaxdzT ! max of dT/dz
REAL(wp), DIMENSION(jpi,jpj) :: zdelr ! delta rho equivalent to deltaT = 0.2
INTEGER :: ji, jj, jk ! dummy loop arguments
REAL(wp) :: zrho3 = 0.03_wp ! density criterion for mixed layer depth
REAL(wp) :: zrho1 = 0.01_wp ! density criterion for mixed layer depth
REAL(wp) :: ztem2 = 0.2_wp ! temperature criterion for mixed layer depth
REAL(wp) :: zztmp, zzdep ! temporary scalars inside do loop
REAL(wp) :: zu, zv, zw, zut, zvt ! temporary workspace
REAL(wp), DIMENSION(A2D(0)) :: zabs2 ! MLD: abs( tn - tn(10m) ) = ztem2
REAL(wp), DIMENSION(A2D(0)) :: ztm2 ! Top of thermocline: tn = tn(10m) - ztem2
REAL(wp), DIMENSION(A2D(0)) :: zrho10_3 ! MLD: rho = rho10m + zrho3
REAL(wp), DIMENSION(A2D(0)) :: zpycn ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
REAL(wp), DIMENSION(A2D(0)) :: ztinv ! max of temperature inversion
REAL(wp), DIMENSION(A2D(0)) :: zdepinv ! depth of temperature inversion
REAL(wp), DIMENSION(A2D(0)) :: zrho0_3 ! MLD rho = rho(surf) = 0.03
REAL(wp), DIMENSION(A2D(0)) :: zrho0_1 ! MLD rho = rho(surf) = 0.01
REAL(wp), DIMENSION(A2D(0)) :: zmaxdzT ! max of dT/dz
REAL(wp), DIMENSION(A2D(0)) :: zdelr ! delta rho equivalent to deltaT = 0.2
!!----------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('dia_hth')
......@@ -131,7 +131,7 @@ CONTAINS
IF( iom_use( 'mlddzt' ) ) zmaxdzT(:,:) = 0._wp
IF( iom_use( 'mlddzt' ) .OR. iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) &
& .OR. iom_use( 'mldr10_3' ) .OR. iom_use( 'pycndep' ) ) THEN
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)
hth (ji,jj) = zztmp
zabs2 (ji,jj) = zztmp
......@@ -142,7 +142,7 @@ CONTAINS
ENDIF
IF( iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) ) THEN
IF( nla10 > 1 ) THEN
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)
zrho0_3(ji,jj) = zztmp
zrho0_1(ji,jj) = zztmp
......@@ -157,7 +157,7 @@ CONTAINS
! MLD: rho = rho(1) + zrho3 !
! MLD: rho = rho(1) + zrho1 !
! ------------------------------------------------------------- !
DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) ! loop from bottom to 2
DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! loop from bottom to 2
!
zzdep = gdepw(ji,jj,jk,Kmm)
zztmp = ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) ) &
......@@ -189,7 +189,7 @@ CONTAINS
!
! Preliminary computation
! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
IF( tmask(ji,jj,nla10) == 1. ) THEN
zu = 1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80 * ts(ji,jj,nla10,jp_sal,Kmm) &
& - 0.0745 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm) &
......@@ -213,7 +213,7 @@ CONTAINS
! temperature inversion: max( 0, max of tn - tn(10m) ) !
! depth of temperature inversion !
! ------------------------------------------------------------- !
DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 ) ! loop from bottom to nlb10
DO_3DS( 0, 0, 0, 0, jpkm1, nlb10, -1 ) ! loop from bottom to nlb10
!
zzdep = gdepw(ji,jj,jk,Kmm) * tmask(ji,jj,1)
!
......@@ -305,13 +305,16 @@ CONTAINS
!
INTEGER :: ji, jj, jk, iid
REAL(wp) :: zztmp, zzdep
INTEGER, DIMENSION(jpi,jpj) :: iktem
INTEGER, DIMENSION(A2D(0)) :: iktem
! --------------------------------------- !
! search deepest level above ptem !
! --------------------------------------- !
iktem(:,:) = 1
DO_3D( 1, 1, 1, 1, 1, jpkm1 ) ! beware temperature is not always decreasing with depth => loop from top to bottom
DO_2D( 0, 0, 0, 0 )
iktem(ji,jj) = 1
END_2D
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! beware temperature is not always decreasing with depth => loop from top to bottom
zztmp = ts(ji,jj,jk,jp_tem,Kmm)
IF( zztmp >= ptem ) iktem(ji,jj) = jk
END_3D
......@@ -319,7 +322,7 @@ CONTAINS
! ------------------------------- !
! Depth of ptem isotherm !
! ------------------------------- !
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
!
zzdep = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ! depth of the ocean bottom
!
......@@ -346,18 +349,29 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phtc
!
INTEGER :: ji, jj, jk, ik
REAL(wp), DIMENSION(jpi,jpj) :: zthick
INTEGER , DIMENSION(jpi,jpj) :: ilevel
REAL(wp), DIMENSION(A2D(0)) :: zthick
INTEGER , DIMENSION(A2D(0)) :: ilevel
! surface boundary condition
IF( .NOT. ln_linssh ) THEN ; zthick(:,:) = 0._wp ; phtc(:,:) = 0._wp
ELSE ; zthick(:,:) = ssh(:,:,Kmm) ; phtc(:,:) = pt(:,:,1) * ssh(:,:,Kmm) * tmask(:,:,1)
IF( .NOT. ln_linssh ) THEN
DO_2D( 0, 0, 0, 0 )
zthick(ji,jj) = 0._wp
phtc (ji,jj) = 0._wp
END_2D
ELSE
DO_2D( 0, 0, 0, 0 )
zthick(ji,jj) = ssh(ji,jj,Kmm)
phtc (ji,jj) = pt(ji,jj,1) * ssh(ji,jj,Kmm) * tmask(ji,jj,1)
END_2D
ENDIF
!
ilevel(:,:) = 1
DO_3D( 1, 1, 1, 1, 1, jpkm1 )
DO_2D( 0, 0, 0, 0 )
ilevel(ji,jj) = 1
END_2D
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
IF( ( gdepw(ji,jj,jk+1,Kmm) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN
ilevel(ji,jj) = jk+1
zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm)
......@@ -365,7 +379,7 @@ CONTAINS
ENDIF
END_3D
!
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
ik = ilevel(ji,jj)
IF( tmask(ji,jj,ik) == 1 ) THEN
zthick(ji,jj) = MIN ( gdepw(ji,jj,ik+1,Kmm), pdep ) - zthick(ji,jj) ! remaining thickness to reach dephw pdep
......
src/OCE/DIA/diamlr.F90
View file @ 78481193
......@@ -6,7 +6,7 @@ MODULE diamlr
!! History : 4.0 ! 2019 (S. Mueller) Original code
!!----------------------------------------------------------------------
USE par_oce , ONLY : wp, jpi, jpj
USE par_oce , ONLY : wp, jpi, jpj, ntsi, ntei, ntsj, ntej
USE phycst , ONLY : rpi
USE dom_oce , ONLY : adatrj
USE tide_mod
......@@ -407,8 +407,9 @@ CONTAINS
!! ** Purpose : update time used in multiple-linear-regression analysis
!!
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj) :: zadatrj2d
REAL(wp), DIMENSION(A2D(0)) :: zadatrj2d
!!----------------------------------------------------------------------
INTEGER :: ji, jj
IF( ln_timing ) CALL timing_start('dia_mlr')
......@@ -417,7 +418,9 @@ CONTAINS
!
! A 2-dimensional field of constant value is sent, and subsequently used directly
! or transformed to a scalar or a constant 3-dimensional field as required.
zadatrj2d(:,:) = adatrj*86400.0_wp
DO_2D( 0, 0, 0, 0 )
zadatrj2d(ji,jj) = adatrj*86400.0_wp
END_2D
IF ( iom_use('diamlr_time') ) CALL iom_put('diamlr_time', zadatrj2d)
!
IF( ln_timing ) CALL timing_stop('dia_mlr')
......
src/OCE/DIA/diaptr.F90
View file @ 78481193
......@@ -78,7 +78,7 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! ocean time-step index
INTEGER , INTENT(in) :: Kmm ! time level index
REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT(in), OPTIONAL :: pvtr ! j-effective transport
REAL(wp), DIMENSION(A2D(0),jpk) , INTENT(in), OPTIONAL :: pvtr ! j-effective transport
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dia_ptr')
......@@ -110,13 +110,13 @@ CONTAINS
!!----------------------------------------------------------------------
!! ** Purpose : Calculate diagnostics and send to XIOS
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! ocean time-step index
INTEGER , INTENT(in) :: Kmm ! time level index
REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT(in), OPTIONAL :: pvtr ! j-effective transport
INTEGER , INTENT(in) :: kt ! ocean time-step index
INTEGER , INTENT(in) :: Kmm ! time level index
REAL(wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL :: pvtr ! j-effective transport (used only by PRESENT)
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace
REAL(wp), DIMENSION(jpj) :: zvsum, ztsum, zssum ! 1D workspace
REAL(wp), DIMENSION(Nis0:Nie0,Njs0:Nje0) :: z2d ! 2D workspace
REAL(wp), DIMENSION( Njs0:Nje0) :: zvsum, ztsum, zssum ! 1D workspace
!
!overturning calculation
REAL(wp), DIMENSION(:,:,: ), ALLOCATABLE :: sjk, r1_sjk, v_msf ! i-mean i-k-surface and its inverse
......@@ -126,19 +126,19 @@ CONTAINS
REAL(wp), DIMENSION(:,:,: ), ALLOCATABLE :: z3dtr
!!----------------------------------------------------------------------
!
ALLOCATE( z3dtr(jpi,jpj,nbasin) )
ALLOCATE( z3dtr(Nis0:Nie0,Njs0:Nje0,nbasin) )
IF( PRESENT( pvtr ) ) THEN
IF( iom_use( 'zomsf' ) ) THEN ! effective MSF
ALLOCATE( z4d1(jpi,jpj,jpk,nbasin) )
ALLOCATE( z4d1(Nis0:Nie0,Njs0:Nje0,jpk,nbasin) )
!
DO jn = 1, nbasin ! by sub-basins
z4d1(1,:,:,jn) = pvtr_int(:,:,jp_vtr,jn) ! zonal cumulative effective transport excluding closed seas
z4d1(Nis0,:,:,jn) = pvtr_int(:,:,jp_vtr,jn) ! zonal cumulative effective transport excluding closed seas
DO jk = jpkm1, 1, -1
z4d1(1,:,jk,jn) = z4d1(1,:,jk+1,jn) - z4d1(1,:,jk,jn) ! effective j-Stream-Function (MSF)
z4d1(Nis0,:,jk,jn) = z4d1(Nis0,:,jk+1,jn) - z4d1(Nis0,:,jk,jn) ! effective j-Stream-Function (MSF)
END DO
DO ji = 2, jpi
z4d1(ji,:,:,jn) = z4d1(1,:,:,jn)
DO ji = Nis0+1, Nie0
z4d1(ji,:,:,jn) = z4d1(Nis0,:,:,jn)
ENDDO
END DO
CALL iom_put( 'zomsf', z4d1 * rc_sv )
......@@ -146,8 +146,8 @@ CONTAINS
DEALLOCATE( z4d1 )
ENDIF
IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN
ALLOCATE( sjk(jpj,jpk,nbasin), r1_sjk(jpj,jpk,nbasin), v_msf(jpj,jpk,nbasin), &
& zt_jk(jpj,jpk,nbasin), zs_jk(jpj,jpk,nbasin) )
ALLOCATE( sjk( Njs0:Nje0,jpk,nbasin), r1_sjk(Njs0:Nje0,jpk,nbasin), v_msf(Njs0:Nje0,jpk,nbasin), &
& zt_jk(Njs0:Nje0,jpk,nbasin), zs_jk( Njs0:Nje0,jpk,nbasin) )
!
DO jn = 1, nbasin
sjk(:,:,jn) = pvtr_int(:,:,jp_msk,jn)
......@@ -162,16 +162,16 @@ CONTAINS
!
ENDDO
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_ove(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_ove(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sophtove', z3dtr )
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_ove(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_ove(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sopstove', z3dtr )
......@@ -181,7 +181,7 @@ CONTAINS
IF( iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN
! Calculate barotropic heat and salt transport here
ALLOCATE( sjk(jpj,1,nbasin), r1_sjk(jpj,1,nbasin) )
ALLOCATE( sjk(A1Dj(0),1,nbasin), r1_sjk(A1Dj(0),1,nbasin) )
!
DO jn = 1, nbasin
sjk(:,1,jn) = SUM( pvtr_int(:,:,jp_msk,jn), 2 )
......@@ -196,16 +196,16 @@ CONTAINS
!
ENDDO
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_btr(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_btr(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sophtbtr', z3dtr )
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_btr(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_btr(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sopstbtr', z3dtr )
......@@ -218,28 +218,28 @@ CONTAINS
pvtr_int(:,:,:,:) = 0._wp
ELSE
IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' ) ) THEN ! i-mean i-k-surface
ALLOCATE( z4d1(jpi,jpj,jpk,nbasin), z4d2(jpi,jpj,jpk,nbasin) )
ALLOCATE( z4d1(Nis0:Nie0,Njs0:Nje0,jpk,nbasin), z4d2(Nis0:Nie0,Njs0:Nje0,jpk,nbasin) )
!
DO jn = 1, nbasin
z4d1(1,:,:,jn) = pzon_int(:,:,jp_msk,jn)
DO ji = 2, jpi
z4d1(ji,:,:,jn) = z4d1(1,:,:,jn)
z4d1(Nis0,:,:,jn) = pzon_int(:,:,jp_msk,jn)
DO ji = Nis0+1, Nie0
z4d1(ji,:,:,jn) = z4d1(Nis0,:,:,jn)
ENDDO
ENDDO
CALL iom_put( 'zosrf', z4d1 )
!
DO jn = 1, nbasin
z4d2(1,:,:,jn) = pzon_int(:,:,jp_tem,jn) / MAX( z4d1(1,:,:,jn), 10.e-15 )
DO ji = 2, jpi
z4d2(ji,:,:,jn) = z4d2(1,:,:,jn)
z4d2(Nis0,:,:,jn) = pzon_int(:,:,jp_tem,jn) / MAX( z4d1(Nis0,:,:,jn), 10.e-15 )
DO ji = Nis0+1, Nie0
z4d2(ji,:,:,jn) = z4d2(Nis0,:,:,jn)
ENDDO
ENDDO
CALL iom_put( 'zotem', z4d2 )
!
DO jn = 1, nbasin
z4d2(1,:,:,jn) = pzon_int(:,:,jp_sal,jn) / MAX( z4d1(1,:,:,jn), 10.e-15 )
DO ji = 2, jpi
z4d2(ji,:,:,jn) = z4d2(1,:,:,jn)
z4d2(Nis0,:,:,jn) = pzon_int(:,:,jp_sal,jn) / MAX( z4d1(Nis0,:,:,jn), 10.e-15 )
DO ji = Nis0+1, Nie0
z4d2(ji,:,:,jn) = z4d2(Nis0,:,:,jn)
ENDDO
ENDDO
CALL iom_put( 'zosal', z4d2 )
......@@ -251,16 +251,16 @@ CONTAINS
IF( iom_use( 'sophtadv' ) .OR. iom_use( 'sopstadv' ) ) THEN
!
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_adv(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_adv(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sophtadv', z3dtr )
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_adv(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_adv(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sopstadv', z3dtr )
......@@ -269,16 +269,16 @@ CONTAINS
IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) THEN
!
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_ldf(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_ldf(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sophtldf', z3dtr )
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_ldf(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_ldf(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sopstldf', z3dtr )
......@@ -287,16 +287,16 @@ CONTAINS
IF( iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) ) THEN
!
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_eiv(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_eiv(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sophteiv', z3dtr )
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_eiv(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_eiv(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sopsteiv', z3dtr )
......@@ -304,16 +304,16 @@ CONTAINS
!
IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_vtr(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_vtr(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sophtvtr', z3dtr )
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_vtr(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_vtr(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sopstvtr', z3dtr )
......@@ -349,8 +349,8 @@ CONTAINS
!! ** Action : pvtr_int - terms for volume streamfunction, heat/salt transport barotropic/overturning terms
!! pzon_int - terms for i mean temperature/salinity
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: Kmm ! time level index
REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in), OPTIONAL :: pvtr ! j-effective transport
INTEGER , INTENT(in) :: Kmm ! time level index
REAL(wp), DIMENSION(A2D(0),jpk), INTENT(in), OPTIONAL :: pvtr ! j-effective transport
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zmask ! 3D workspace
REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zts ! 4D workspace
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: sjk, v_msf ! Zonal sum: i-k surface area, j-effective transport
......@@ -362,7 +362,7 @@ CONTAINS
IF( PRESENT( pvtr ) ) THEN
! i sum of effective j transport excluding closed seas
IF( iom_use( 'zomsf' ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN
ALLOCATE( v_msf(A1Dj(nn_hls),jpk,nbasin) )
ALLOCATE( v_msf(A1Dj(0),jpk,nbasin) )
DO jn = 1, nbasin
v_msf(:,:,jn) = ptr_sjk( pvtr(:,:,:), btmsk34(:,:,jn) )
......@@ -374,16 +374,16 @@ CONTAINS
ENDIF
! i sum of j surface area, j surface area - temperature/salinity product on V grid
IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR. &
IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR. &
& iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN
ALLOCATE( zmask(A2D(nn_hls),jpk), zts(A2D(nn_hls),jpk,jpts), &
& sjk(A1Dj(nn_hls),jpk,nbasin), &
& zt_jk(A1Dj(nn_hls),jpk,nbasin), zs_jk(A1Dj(nn_hls),jpk,nbasin) )
ALLOCATE( zmask( A2D(0),jpk ), zts( A2D(0),jpk,jpts ), &
& sjk( A1Dj(0),jpk,nbasin), &
& zt_jk(A1Dj(0),jpk,nbasin), zs_jk(A1Dj(0),jpk,nbasin) )
zmask(:,:,:) = 0._wp
zts(:,:,:,:) = 0._wp
DO_3D( 1, 1, 1, 0, 1, jpkm1 )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
zmask(ji,jj,jk) = vmask(ji,jj,jk) * zvfc
zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc !Tracers averaged onto V grid
......@@ -405,14 +405,14 @@ CONTAINS
ELSE
! i sum of j surface area - temperature/salinity product on T grid
IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' ) ) THEN
ALLOCATE( zmask(A2D(nn_hls),jpk), zts(A2D(nn_hls),jpk,jpts), &
& sjk(A1Dj(nn_hls),jpk,nbasin), &
& zt_jk(A1Dj(nn_hls),jpk,nbasin), zs_jk(A1Dj(nn_hls),jpk,nbasin) )
ALLOCATE( zmask( A2D(0),jpk ), zts( A2D(0),jpk,jpts ), &
& sjk( A1Dj(0),jpk,nbasin), &
& zt_jk(A1Dj(0),jpk,nbasin), zs_jk(A1Dj(0),jpk,nbasin) )
zmask(:,:,:) = 0._wp
zts(:,:,:,:) = 0._wp
DO_3D( 1, 1, 1, 1, 1, jpkm1 )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zsfc = e1t(ji,jj) * e3t(ji,jj,jk,Kmm)
zmask(ji,jj,jk) = tmask(ji,jj,jk) * zsfc
zts(ji,jj,jk,jp_tem) = ts(ji,jj,jk,jp_tem,Kmm) * zsfc
......@@ -434,11 +434,12 @@ CONTAINS
! i-k sum of j surface area - temperature/salinity product on V grid
IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN
! TODO: Can be A2D(0) once all dia_ptr_hst calls have arguments with consistent declarations
ALLOCATE( zts(A2D(nn_hls),jpk,jpts) )
zts(:,:,:,:) = 0._wp
DO_3D( 1, 1, 1, 0, 1, jpkm1 )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc !Tracers averaged onto V grid
zts(ji,jj,jk,jp_sal) = (ts(ji,jj,jk,jp_sal,Kmm)+ts(ji,jj+1,jk,jp_sal,Kmm)) * 0.5 * zvfc
......@@ -538,11 +539,12 @@ CONTAINS
!! Wrapper for heat and salt transport calculations to calculate them for each basin
!! Called from all advection and/or diffusion routines
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ktra ! tracer index
CHARACTER(len=3) , INTENT(in) :: cptr ! transport type 'adv'/'ldf'/'eiv'
REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT(in) :: pvflx ! 3D input array of advection/diffusion
REAL(wp), DIMENSION(A1Dj(nn_hls),nbasin) :: zsj !
INTEGER :: jn !
INTEGER , INTENT(in) :: ktra ! tracer index
CHARACTER(len=3) , INTENT(in) :: cptr ! transport type 'adv'/'ldf'/'eiv'
! TODO: Can be A2D(0) once all dia_ptr_hst calls have arguments with consistent declarations
REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pvflx ! 3D input array of advection/diffusion
REAL(wp), DIMENSION(A1Dj(0),nbasin) :: zsj !
INTEGER :: jn !
DO jn = 1, nbasin
zsj(:,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
......@@ -576,13 +578,13 @@ CONTAINS
!!
!! ** Action : phstr
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpj,nbasin) , INTENT(inout) :: phstr !
REAL(wp), DIMENSION(A1Dj(nn_hls),nbasin), INTENT(in) :: pva !
REAL(wp), DIMENSION(Njs0:Nje0,nbasin), INTENT(inout) :: phstr !
REAL(wp), DIMENSION(A1Dj(0) ,nbasin), INTENT(in ) :: pva !
INTEGER :: jj
#if ! defined key_mpi_off
INTEGER, DIMENSION(1) :: ish1d
INTEGER, DIMENSION(2) :: ish2d
REAL(wp), DIMENSION(jpj*nbasin) :: zwork
INTEGER, DIMENSION(1) :: ish1d
INTEGER, DIMENSION(2) :: ish2d
REAL(wp), DIMENSION(:), ALLOCATABLE :: zwork
#endif
DO jj = ntsj, ntej
......@@ -591,11 +593,13 @@ CONTAINS
#if ! defined key_mpi_off
IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN
ish1d(1) = jpj*nbasin
ish2d(1) = jpj ; ish2d(2) = nbasin
ALLOCATE( zwork(Nj_0*nbasin) )
ish1d(1) = Nj_0*nbasin
ish2d(1) = Nj_0 ; ish2d(2) = nbasin
zwork(:) = RESHAPE( phstr(:,:), ish1d )
CALL mpp_sum( 'diaptr', zwork, ish1d(1), ncomm_znl )
phstr(:,:) = RESHAPE( zwork, ish2d )
DEALLOCATE( zwork )
ENDIF
#endif
END SUBROUTINE ptr_sum_2d
......@@ -612,13 +616,13 @@ CONTAINS
!!
!! ** Action : phstr
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpj,jpk,nbasin) , INTENT(inout) :: phstr !
REAL(wp), DIMENSION(A1Dj(nn_hls),jpk,nbasin), INTENT(in) :: pva !
INTEGER :: jj, jk
REAL(wp), DIMENSION(Njs0:Nje0,jpk,nbasin), INTENT(inout) :: phstr !
REAL(wp), DIMENSION(A1Dj(0) ,jpk,nbasin), INTENT(in ) :: pva !
INTEGER :: jj, jk
#if ! defined key_mpi_off
INTEGER, DIMENSION(1) :: ish1d
INTEGER, DIMENSION(3) :: ish3d
REAL(wp), DIMENSION(jpj*jpk*nbasin) :: zwork
INTEGER, DIMENSION(1) :: ish1d
INTEGER, DIMENSION(3) :: ish3d
REAL(wp), DIMENSION(:), ALLOCATABLE :: zwork
#endif
DO jk = 1, jpk
......@@ -629,11 +633,13 @@ CONTAINS
#if ! defined key_mpi_off
IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN
ish1d(1) = jpj*jpk*nbasin
ish3d(1) = jpj ; ish3d(2) = jpk ; ish3d(3) = nbasin
ALLOCATE( zwork(Nj_0*jpk*nbasin) )
ish1d(1) = Nj_0*jpk*nbasin
ish3d(1) = Nj_0 ; ish3d(2) = jpk ; ish3d(3) = nbasin
zwork(:) = RESHAPE( phstr(:,:,:), ish1d )
CALL mpp_sum( 'diaptr', zwork, ish1d(1), ncomm_znl )
phstr(:,:,:) = RESHAPE( zwork, ish3d )
DEALLOCATE( zwork )
ENDIF
#endif
END SUBROUTINE ptr_sum_3d
......@@ -651,13 +657,13 @@ CONTAINS
! nbasin has been initialized in iom_init to define the axis "basin"
!
IF( .NOT. ALLOCATED( btmsk ) ) THEN
ALLOCATE( btmsk(jpi,jpj,nbasin) , btmsk34(jpi,jpj,nbasin), &
& hstr_adv(jpj,jpts,nbasin), hstr_eiv(jpj,jpts,nbasin), &
& hstr_ove(jpj,jpts,nbasin), hstr_btr(jpj,jpts,nbasin), &
& hstr_ldf(jpj,jpts,nbasin), hstr_vtr(jpj,jpts,nbasin), STAT=ierr(1) )
ALLOCATE( btmsk(jpi,jpj,nbasin) , btmsk34(jpi,jpj,nbasin), &
& hstr_adv(Njs0:Nje0,jpts,nbasin), hstr_eiv(Njs0:Nje0,jpts,nbasin), &
& hstr_ove(Njs0:Nje0,jpts,nbasin), hstr_btr(Njs0:Nje0,jpts,nbasin), &
& hstr_ldf(Njs0:Nje0,jpts,nbasin), hstr_vtr(Njs0:Nje0,jpts,nbasin), STAT=ierr(1) )
!
ALLOCATE( pvtr_int(jpj,jpk,jpts+2,nbasin), &
& pzon_int(jpj,jpk,jpts+1,nbasin), STAT=ierr(2) )
ALLOCATE( pvtr_int(Njs0:Nje0,jpk,jpts+2,nbasin), &
& pzon_int(Njs0:Nje0,jpk,jpts+1,nbasin), STAT=ierr(2) )
!
dia_ptr_alloc = MAXVAL( ierr )
CALL mpp_sum( 'diaptr', dia_ptr_alloc )
......@@ -677,11 +683,12 @@ CONTAINS
!!
!! ** Action : - p_fval: i-k-mean poleward flux of pvflx
!!----------------------------------------------------------------------
! TODO: Can be A2D(0) once all dia_ptr_hst calls have arguments with consistent declarations
REAL(wp), INTENT(in), DIMENSION(A2D(nn_hls),jpk) :: pvflx ! mask flux array at V-point
REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: pmsk ! Optional 2D basin mask
REAL(wp), INTENT(in), DIMENSION(jpi,jpj ) :: pmsk ! Optional 2D basin mask
!
INTEGER :: ji, jj, jk ! dummy loop arguments
REAL(wp), DIMENSION(A1Dj(nn_hls)) :: p_fval ! function value
INTEGER :: ji, jj, jk ! dummy loop arguments
REAL(wp), DIMENSION(A1Dj(0)) :: p_fval ! function value
!!--------------------------------------------------------------------
!
p_fval(:) = 0._wp
......@@ -702,11 +709,12 @@ CONTAINS
!!
!! ** Action : - p_fval: i-k-mean poleward flux of pvflx
!!----------------------------------------------------------------------
REAL(wp) , INTENT(in), DIMENSION(A2D(nn_hls)) :: pvflx ! mask flux array at V-point
REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) :: pmsk ! Optional 2D basin mask
! TODO: Can be A2D(0) once all dia_ptr_hst calls have arguments with consistent declarations
REAL(wp) , INTENT(in), DIMENSION(A2D(nn_hls)) :: pvflx ! mask flux array at V-point
REAL(wp) , INTENT(in), DIMENSION(jpi,jpj ) :: pmsk ! Optional 2D basin mask
!
INTEGER :: ji,jj ! dummy loop arguments
REAL(wp), DIMENSION(A1Dj(nn_hls)) :: p_fval ! function value
INTEGER :: ji, jj ! dummy loop arguments
REAL(wp), DIMENSION(A1Dj(0)) :: p_fval ! function value
!!--------------------------------------------------------------------
!
p_fval(:) = 0._wp
......@@ -725,14 +733,13 @@ CONTAINS
!!
!! ** Action : - p_fval: j-cumulated sum of pva
!!----------------------------------------------------------------------
REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) :: pva ! mask flux array at V-point
REAL(wp) , INTENT(in), DIMENSION(A2D(0)) :: pva ! mask flux array at V-point
!
INTEGER :: ji,jj,jc ! dummy loop arguments
INTEGER :: ijpj ! ???
REAL(wp), DIMENSION(jpi,jpj) :: p_fval ! function value
INTEGER :: ji,jj,jc ! dummy loop arguments
INTEGER :: ijpj ! ???
REAL(wp), DIMENSION(A2D(0)) :: p_fval ! function value
!!--------------------------------------------------------------------
!
ijpj = jpj ! ???
p_fval(:,:) = 0._wp
DO jc = 1, jpnj ! looping over all processors in j axis
DO_2D( 0, 0, 0, 0 )
......@@ -756,11 +763,11 @@ CONTAINS
!!----------------------------------------------------------------------
!!
IMPLICIT none
REAL(wp) , INTENT(in), DIMENSION(A2D(nn_hls),jpk) :: pta ! mask flux array at V-point
REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) :: pmsk ! Optional 2D basin mask
REAL(wp) , INTENT(in), DIMENSION(A2D(0) ,jpk) :: pta ! mask flux array at V-point
REAL(wp) , INTENT(in), DIMENSION(jpi,jpj ) :: pmsk ! Optional 2D basin mask
!!
INTEGER :: ji, jj, jk ! dummy loop arguments
REAL(wp), DIMENSION(A1Dj(nn_hls),jpk) :: p_fval ! return function value
REAL(wp), DIMENSION(A1Dj(0),jpk) :: p_fval ! return function value
!!--------------------------------------------------------------------
!
p_fval(:,:) = 0._wp
......
src/OCE/DIA/diawri.F90
View file @ 78481193
......@@ -135,8 +135,8 @@ CONTAINS
ENDIF
! initialize arrays
z2d(:,:) = 0._wp
z3d(:,:,:) = 0._wp
z2d(A2D(0)) = 0._wp
z3d(A2D(0),:) = 0._wp
! Output of initial vertical scale factor
CALL iom_put("e3t_0", e3t_0(:,:,:) )
......@@ -868,7 +868,11 @@ CONTAINS
CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "J/m2" , & ! htc3
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
#endif
CALL histdef( nid_T, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
& jpi, jpj, nh_T, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
CALL histdef( nid_T, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
& jpi, jpj, nh_T, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
!
CALL histend( nid_T, snc4chunks=snc4set )
! !!! nid_U : 3D
......@@ -878,10 +882,7 @@ CONTAINS
CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
ENDIF
! !!! nid_U : 2D
CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
!
CALL histend( nid_U, snc4chunks=snc4set )
! !!! nid_V : 3D
......@@ -891,10 +892,7 @@ CONTAINS
CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
ENDIF
! !!! nid_V : 2D
CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
!
CALL histend( nid_V, snc4chunks=snc4set )
! !!! nid_W : 3D
......@@ -1066,12 +1064,12 @@ CONTAINS
CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
#endif
CALL histwrite( nid_T, "sozotaux", it, utau , ndim_hT, ndex_hT ) ! i-wind stress
CALL histwrite( nid_T, "sometauy", it, vtau , ndim_hT, ndex_hT ) ! j-wind stress
CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm) , ndim_U , ndex_U ) ! i-current
CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm) , ndim_V , ndex_V ) ! j-current
CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
IF( ln_zad_Aimp ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
......
src/OCE/DYN/dynatf.F90
View file @ 78481193
......@@ -331,7 +331,7 @@ CONTAINS
ALLOCATE(zutau(jpi,jpj))
DO_2D( 0, 0, 0, 0 )
jk = miku(ji,jj)
zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * puu(ji,jj,jk,Kaa)
zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( puu(ji-1,jj,jk,Kaa) + puu(ji,jj,jk,Kaa) )
END_2D
CALL iom_put( "utau", zutau(:,:) )
DEALLOCATE(zutau)
......@@ -345,7 +345,7 @@ CONTAINS
ALLOCATE(zvtau(jpi,jpj))
DO_2D( 0, 0, 0, 0 )
jk = mikv(ji,jj)
zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * pvv(ji,jj,jk,Kaa)
zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( pvv(ji,jj-1,jk,Kaa) + pvv(ji,jj,jk,Kaa) )
END_2D
CALL iom_put( "vtau", zvtau(:,:) )
DEALLOCATE(zvtau)
......
src/OCE/DYN/dynatf_qco.F90
View file @ 78481193
......@@ -248,7 +248,7 @@ CONTAINS
ALLOCATE(zutau(jpi,jpj))
DO_2D( 0, 0, 0, 0 )
jk = miku(ji,jj)
zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * puu(ji,jj,jk,Kaa)
zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( puu(ji-1,jj,jk,Kaa) + puu(ji,jj,jk,Kaa) )
END_2D
CALL iom_put( "utau", zutau(:,:) )
DEALLOCATE(zutau)
......@@ -262,7 +262,7 @@ CONTAINS
ALLOCATE(zvtau(jpi,jpj))
DO_2D( 0, 0, 0, 0 )
jk = mikv(ji,jj)
zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * pvv(ji,jj,jk,Kaa)
zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( pvv(ji,jj-1,jk,Kaa) + pvv(ji,jj,jk,Kaa) )
END_2D
CALL iom_put( "vtau", zvtau(:,:) )
DEALLOCATE(zvtau)
......
src/OCE/DYN/dynspg_ts.F90
View file @ 78481193
......@@ -334,14 +334,14 @@ CONTAINS
! ! ------------------ !
IF( ln_bt_fw ) THEN
DO_2D( 0, 0, 0, 0 )
zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm)
zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm)
zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utauU(ji,jj) * r1_hu(ji,jj,Kmm)
zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtauV(ji,jj) * r1_hv(ji,jj,Kmm)
END_2D
ELSE
zztmp = r1_rho0 * r1_2
DO_2D( 0, 0, 0, 0 )
zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm)
zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm)
zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utauU(ji,jj) ) * r1_hu(ji,jj,Kmm)
zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtauV(ji,jj) ) * r1_hv(ji,jj,Kmm)
END_2D
ENDIF
!
......
src/OCE/DYN/dynzdf.F90
View file @ 78481193
......@@ -267,10 +267,10 @@ CONTAINS
DO_2D( 0, 0, 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
#if defined key_RK3
! ! RK3: use only utau (not utau_b)
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * utau(ji,jj) &
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * utauU(ji,jj) &
& / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
#else
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utau(ji,jj) ) &
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utauU(ji,jj) ) &
& / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
#endif
END_2D
......@@ -397,10 +397,10 @@ CONTAINS
DO_2D( 0, 0, 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
#if defined key_RK3
! ! RK3: use only vtau (not vtau_b)
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * vtau(ji,jj) &
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * vtauV(ji,jj) &
& / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
#else
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2*( vtau_b(ji,jj) + vtau(ji,jj) ) &
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2 * ( vtau_b(ji,jj) + vtauV(ji,jj) ) &
& / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
#endif
END_2D
......
src/OCE/IOM/iom.F90
View file @ 78481193
......@@ -1327,14 +1327,21 @@ CONTAINS
IF( irankpv == 1 ) ishape(1:1) = SHAPE(pv_r1d)
IF( irankpv == 2 ) ishape(1:2) = SHAPE(pv_r2d)
IF( irankpv == 3 ) ishape(1:3) = SHAPE(pv_r3d)
ix1 = 1 ; ix2 = icnt(1) ; iy1 = 1 ; iy2 = icnt(2) ! index of the array to be read
ctmp1 = 'd'
ELSE
IF( irankpv == 2 ) THEN
ishape(1:2) = SHAPE(pv_r2d(Nis0:Nie0,Njs0:Nje0 )) ; ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0)'
ENDIF
IF( irankpv == 3 ) THEN
ishape(1:3) = SHAPE(pv_r3d(Nis0:Nie0,Njs0:Nje0,:)) ; ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0,:)'
IF( irankpv == 2 ) ishape(1:2) = SHAPE(pv_r2d)
IF( irankpv == 3 ) ishape(1:3) = SHAPE(pv_r3d)
IF( ishape(1) == jpi .AND. ishape(2) == jpj ) THEN
ishape(1:2) = (/ Ni_0, Nj_0 /)
ix1 = Nis0 ; ix2 = Nie0 ; iy1 = Njs0 ; iy2 = Nje0 ! index of the array to be read
ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0'
ELSE
ix1 = 1 ; ix2 = icnt(1) ; iy1 = 1 ; iy2 = icnt(2) ! index of the array to be read
ctmp1 = 'd(:,:'
ENDIF
IF( irankpv == 3 ) ctmp1 = TRIM(ctmp1)//',:'
ctmp1 = TRIM(ctmp1)//')'
ENDIF
DO jl = 1, irankpv
WRITE( ctmp2, FMT="(', ', i1,'): ', i5,' /= icnt(', i1,'):', i5)" ) jl, ishape(jl), jl, icnt(jl)
......@@ -1347,11 +1354,6 @@ CONTAINS
!-
IF( idvar > 0 .AND. istop == nstop ) THEN ! no additional errors until this point...
!
! find the right index of the array to be read
IF( idom /= jpdom_unknown ) THEN ; ix1 = Nis0 ; ix2 = Nie0 ; iy1 = Njs0 ; iy2 = Nje0
ELSE ; ix1 = 1 ; ix2 = icnt(1) ; iy1 = 1 ; iy2 = icnt(2)
ENDIF
CALL iom_nf90_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d, pv_r3d )
IF( istop == nstop ) THEN ! no additional errors until this point...
......@@ -1362,9 +1364,11 @@ CONTAINS
zsgn = 1._wp
IF( PRESENT(psgn ) ) zsgn = psgn
!--- overlap areas and extra hallows (mpp)
IF( PRESENT(pv_r2d) .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
llok = idom /= jpdom_unknown .AND. cl_type /= 'Z' &
& .AND. ix1 == Nis0 .AND. ix2 == Nie0 .AND. iy1 == Njs0 .AND. iy2 == Nje0
IF( PRESENT(pv_r2d) .AND. llok ) THEN
CALL lbc_lnk( 'iom', pv_r2d, cl_type, zsgn, kfillmode = kfill )
ELSEIF( PRESENT(pv_r3d) .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
ELSEIF( PRESENT(pv_r3d) .AND. llok ) THEN
CALL lbc_lnk( 'iom', pv_r3d, cl_type, zsgn, kfillmode = kfill )
ENDIF
!
......@@ -2336,14 +2340,14 @@ CONTAINS
idb(jn) = -nn_hls ! Tile data offset (halo size)
END DO
! Tile_[ij]begin are defined with respect to the processor data domain, so data_[ij]begin is added
CALL iom_set_domain_attr("grid_"//cdgrd, ntiles=nijtile, &
& tile_ibegin=ntsi_a(1:nijtile) + idb(:) - 1, tile_jbegin=ntsj_a(1:nijtile) + idb(:) - 1, &
& tile_ibegin=ntsi_a(1:nijtile) - nn_hls - 1, tile_jbegin=ntsj_a(1:nijtile) - nn_hls - 1, &
& tile_ni=ini(:), tile_nj=inj(:), &
& tile_data_ibegin=idb(:), tile_data_jbegin=idb(:), &
& tile_data_ni=ini(:) - 2 * idb(:), tile_data_nj=inj(:) - 2 * idb(:))
idb(:) = 0
CALL iom_set_domain_attr("grid_"//cdgrd//"_inner", ntiles=nijtile, &
& tile_ibegin=ntsi_a(1:nijtile) + idb(:) - 1, tile_jbegin=ntsj_a(1:nijtile) + idb(:) - 1, &
& tile_ibegin=ntsi_a(1:nijtile) - nn_hls - 1, tile_jbegin=ntsj_a(1:nijtile) - nn_hls - 1, &
& tile_ni=ini(:), tile_nj=inj(:), &
& tile_data_ibegin=idb(:), tile_data_jbegin=idb(:), &
& tile_data_ni=ini(:) - 2 * idb(:), tile_data_nj=inj(:) - 2 * idb(:))
......@@ -2453,7 +2457,7 @@ CONTAINS
! CALL dom_ngb( -168.53_wp, 65.03_wp, ix, iy, 'T' ) ! i-line that passes through Bering Strait: Reference latitude (used in plots)
CALL dom_ngb( 180.0_wp, 90.0_wp, ix, iy, 'T' ) ! i-line that passes near the North Pole : Reference latitude (used in plots)
CALL iom_set_domain_attr("gznl", ni_glo=Ni0glo, nj_glo=Nj0glo, ibegin=mig0(Nis0)-1, jbegin=mjg0(Njs0)-1, ni=Ni_0, nj=Nj_0)
CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = -nn_hls, data_ni = jpi, data_jbegin = -nn_hls, data_nj = jpj)
CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin=0, data_ni=Ni_0, data_jbegin=0, data_nj=Nj_0)
CALL iom_set_domain_attr("gznl", lonvalue = real(zlon, dp), &
& latvalue = real(RESHAPE(plat(Nis0:Nie0, Njs0:Nje0),(/ Ni_0*Nj_0 /)),dp))
CALL iom_set_zoom_domain_attr("ptr", ibegin=ix-1, jbegin=0, ni=1, nj=Nj0glo)
......
src/OCE/IOM/iom_nf90.F90
View file @ 78481193
......@@ -40,6 +40,8 @@ MODULE iom_nf90
MODULE PROCEDURE iom_nf90_rp0123d_dp
END INTERFACE
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: iom_nf90.F90 14433 2021-02-11 08:06:49Z smasson $
......@@ -544,7 +546,7 @@ CONTAINS
INTEGER :: idvar ! variable id
INTEGER :: jd ! dimension loop counter
INTEGER :: ix1, ix2, iy1, iy2 ! subdomain indexes
INTEGER, DIMENSION(4) :: idimsz ! dimensions size
INTEGER, DIMENSION(3) :: ishape ! dimensions size
INTEGER, DIMENSION(4) :: idimid ! dimensions id
CHARACTER(LEN=256) :: clinfo ! info character
INTEGER :: if90id ! nf90 file identifier
......@@ -627,11 +629,9 @@ CONTAINS
itype = NF90_DOUBLE
ENDIF
IF( PRESENT(pv_r0d) ) THEN
CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, &
& iom_file(kiomid)%nvid(idvar) ), clinfo )
CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, iom_file(kiomid)%nvid(idvar) ), clinfo )
ELSE
CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, idimid(1:idims), &
& iom_file(kiomid)%nvid(idvar) ), clinfo )
CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, idimid(1:idims), iom_file(kiomid)%nvid(idvar) ), clinfo )
ENDIF
lchunk = .false.
IF( snc4set%luse .AND. idims == 4 ) lchunk = .true.
......@@ -673,23 +673,13 @@ CONTAINS
ENDIF
! on what kind of domain must the data be written?
IF( PRESENT(pv_r2d) .OR. PRESENT(pv_r3d) ) THEN
idimsz(1:2) = iom_file(kiomid)%dimsz(1:2,idvar)
IF( idimsz(1) == Ni_0 .AND. idimsz(2) == Nj_0 ) THEN
ix1 = Nis0 ; ix2 = Nie0 ; iy1 = Njs0 ; iy2 = Nje0
ELSEIF( idimsz(1) == jpi .AND. idimsz(2) == jpj ) THEN
ix1 = 1 ; ix2 = jpi ; iy1 = 1 ; iy2 = jpj
ELSEIF( idimsz(1) == jpi .AND. idimsz(2) == jpj ) THEN
ix1 = 1 ; ix2 = jpi ; iy1 = 1 ; iy2 = jpj
ELSE
CALL ctl_stop( 'iom_nf90_rp0123d: should have been an impossible case...' )
ENDIF
! write dimension variables if it is not already done
! =============
! trick: is defined to 0 => dimension variable are defined but not yet written
IF( iom_file(kiomid)%dimsz(1, 4) == 0 ) THEN ! time_counter = 0
CALL iom_nf90_check( NF90_PUT_VAR( if90id, 1, glamt(ix1:ix2, iy1:iy2) ), clinfo )
CALL iom_nf90_check( NF90_PUT_VAR( if90id, 2, gphit(ix1:ix2, iy1:iy2) ), clinfo )
CALL iom_nf90_check( NF90_PUT_VAR( if90id, 1, glamt(A2D(0)) ), clinfo )
CALL iom_nf90_check( NF90_PUT_VAR( if90id, 2, gphit(A2D(0)) ), clinfo )
SELECT CASE (iom_file(kiomid)%comp)
CASE ('OCE')
CALL iom_nf90_check( NF90_PUT_VAR( if90id, 3, gdept_1d ), clinfo )
......@@ -704,6 +694,17 @@ CONTAINS
iom_file(kiomid)%dimsz(1, 4) = 1 ! so we don't enter this IF case any more...
IF(lwp) WRITE(numout,*) TRIM(clinfo)//' write dimension variables done'
ENDIF
IF( PRESENT(pv_r2d) ) ishape(1:2) = SHAPE(pv_r2d)
IF( PRESENT(pv_r3d) ) ishape(1:3) = SHAPE(pv_r3d)
IF( ishape(1) == Ni_0 .AND. ishape(2) == Nj_0 ) THEN
ix1 = 1 ; ix2 = Ni_0 ; iy1 = 1 ; iy2 = Nj_0
ELSEIF( ishape(1) == jpi .AND. ishape(2) == jpj ) THEN
ix1 = Nis0 ; ix2 = Nie0 ; iy1 = Njs0 ; iy2 = Nje0
ELSE
CALL ctl_stop( 'iom_nf90_rp0123d: should have been an impossible case...' )
ENDIF
ENDIF
! write the data
......@@ -712,7 +713,7 @@ CONTAINS
CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r0d ), clinfo )
ELSEIF( PRESENT(pv_r1d) ) THEN
CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r1d(:) ), clinfo )
ELSEIF( PRESENT(pv_r2d) ) THEN
ELSEIF( PRESENT(pv_r2d) ) THEN
CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r2d(ix1:ix2,iy1:iy2) ), clinfo )
ELSEIF( PRESENT(pv_r3d) ) THEN
CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r3d(ix1:ix2,iy1:iy2,:) ), clinfo )
......
src/OCE/IOM/prtctl.F90
View file @ 78481193
......@@ -137,9 +137,14 @@ CONTAINS
INTEGER :: jn, jl, kdir
INTEGER :: iis, iie, jjs, jje
INTEGER :: itra, inum
INTEGER, DIMENSION(4) :: ishape
REAL(2*wp) :: zsum1, zsum2, zvctl1, zvctl2
!!----------------------------------------------------------------------
!
IF( ( ktab2d_1 * ktab3d_1 * ktab4d_1 * ktab2d_2 * ktab3d_2 ) /= 0 ) THEN
CALL ctl_stop( 'prt_ctl is not working with tiles' )
ENDIF
! Arrays, scalars initialization
cl1 = ''
cl2 = ''
......@@ -157,12 +162,19 @@ CONTAINS
! Loop over each sub-domain, i.e. the total number of processors ijsplt
DO jl = 1, SIZE(nall_ictls)
! define shoter names...
iis = MAX( nall_ictls(jl), ntsi )
iie = MIN( nall_ictle(jl), ntei )
jjs = MAX( nall_jctls(jl), ntsj )
jje = MIN( nall_jctle(jl), ntej )
IF( PRESENT(tab2d_1) ) ishape(1:2) = SHAPE(tab2d_1)
IF( PRESENT(tab3d_1) ) ishape(1:3) = SHAPE(tab3d_1)
IF( PRESENT(tab4d_1) ) ishape(1:4) = SHAPE(tab4d_1)
IF( ishape(1) == jpi .AND. ishape(2) == jpj ) THEN
iis = Nis0 ; iie = Nie0 ; jjs = Njs0 ; jje = Nje0
ELSE
iis = 1 ; iie = ishape(1) ; jjs = 1 ; jje = ishape(2)
ENDIF
iis = MAX( nall_ictls(jl), iis )
iie = MIN( nall_ictle(jl), iie )
jjs = MAX( nall_jctls(jl), jjs )
jje = MIN( nall_jctle(jl), jje )
IF( PRESENT(clinfo) ) THEN ; inum = numprt_top(jl)
ELSE ; inum = numprt_oce(jl)
......@@ -188,32 +200,32 @@ CONTAINS
! 2D arrays
IF( PRESENT(tab2d_1) ) THEN
IF( PRESENT(mask1) ) THEN ; zsum1 = SUM( tab2d_1(iis:iie,jjs:jje) * mask1(iis:iie,jjs:jje,1) )
ELSE ; zsum1 = SUM( tab2d_1(iis:iie,jjs:jje) )
IF( PRESENT(mask1) ) THEN ; zsum1 = SUM( tab2d_1(iis:iie,jjs:jje) * mask1(A2D(0),1) )
ELSE ; zsum1 = SUM( tab2d_1(iis:iie,jjs:jje) )
ENDIF
ENDIF
IF( PRESENT(tab2d_2) ) THEN
IF( PRESENT(mask2) ) THEN ; zsum2 = SUM( tab2d_2(iis:iie,jjs:jje) * mask2(iis:iie,jjs:jje,1) )
ELSE ; zsum2 = SUM( tab2d_2(iis:iie,jjs:jje) )
IF( PRESENT(mask2) ) THEN ; zsum2 = SUM( tab2d_2(iis:iie,jjs:jje) * mask2(A2D(0),1) )
ELSE ; zsum2 = SUM( tab2d_2(iis:iie,jjs:jje) )
ENDIF
ENDIF
! 3D arrays
IF( PRESENT(tab3d_1) ) THEN
IF( PRESENT(mask1) ) THEN ; zsum1 = SUM( tab3d_1(iis:iie,jjs:jje,1:kdir) * mask1(iis:iie,jjs:jje,1:kdir) )
ELSE ; zsum1 = SUM( tab3d_1(iis:iie,jjs:jje,1:kdir) )
IF( PRESENT(mask1) ) THEN ; zsum1 = SUM( tab3d_1(iis:iie,jjs:jje,1:kdir) * mask1(A2D(0),1:kdir) )
ELSE ; zsum1 = SUM( tab3d_1(iis:iie,jjs:jje,1:kdir) )
ENDIF
ENDIF
IF( PRESENT(tab3d_2) ) THEN
IF( PRESENT(mask2) ) THEN ; zsum2 = SUM( tab3d_2(iis:iie,jjs:jje,1:kdir) * mask2(iis:iie,jjs:jje,1:kdir) )
ELSE ; zsum2 = SUM( tab3d_2(iis:iie,jjs:jje,1:kdir) )
IF( PRESENT(mask2) ) THEN ; zsum2 = SUM( tab3d_2(iis:iie,jjs:jje,1:kdir) * mask2(A2D(0),1:kdir) )
ELSE ; zsum2 = SUM( tab3d_2(iis:iie,jjs:jje,1:kdir) )
ENDIF
ENDIF
! 4D arrays
IF( PRESENT(tab4d_1) ) THEN
IF( PRESENT(mask1) ) THEN ; zsum1 = SUM( tab4d_1(iis:iie,jjs:jje,1:kdir,jn) * mask1(iis:iie,jjs:jje,1:kdir) )
ELSE ; zsum1 = SUM( tab4d_1(iis:iie,jjs:jje,1:kdir,jn) )
IF( PRESENT(mask1) ) THEN ; zsum1 = SUM( tab4d_1(iis:iie,jjs:jje,1:kdir,jn) * mask1(A2D(0),1:kdir) )
ELSE ; zsum1 = SUM( tab4d_1(iis:iie,jjs:jje,1:kdir,jn) )
ENDIF
ENDIF
......
src/OCE/SBC/cpl_oasis3.F90
View file @ 78481193
......@@ -419,6 +419,7 @@ CONTAINS
IF( .NOT. ll_1st ) THEN
CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,jc), srcv(kid)%clgrid, srcv(kid)%nsgn )
ENDIF
!!clem: mettre T instead of clgrid
ENDDO
!
......
src/OCE/SBC/sbc_ice.F90
View file @ 78481193
......@@ -50,8 +50,8 @@ MODULE sbc_ice
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qcn_ice !: heat conduction flux in the layer below surface [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qtr_ice_top !: solar flux transmitted below the ice surface [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau_ice !: atmos-ice u-stress. VP: I-pt ; EVP: U,V-pts [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtau_ice !: atmos-ice v-stress. VP: I-pt ; EVP: U,V-pts [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau_ice !: atmos-ice u-stress. T-pts [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtau_ice !: atmos-ice v-stress. T-pts [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: emp_ice !: sublimation - precip over sea ice [kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: topmelt !: category topmelt
......
src/OCE/SBC/sbc_oce.F90
View file @ 78481193
......@@ -103,34 +103,36 @@ MODULE sbc_oce
INTEGER , PUBLIC :: ncpl_qsr_freq = 0 !: qsr coupling frequency per days from atmosphere (used by top)
!
!! !! now ! before !!
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau , utau_b !: sea surface i-stress (ocean referential) [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtau , vtau_b !: sea surface j-stress (ocean referential) [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau_icb, vtau_icb !: sea surface (i,j)-stress used by icebergs [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: taum !: module of sea surface stress (at T-point) [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau !: sea surface i-stress (ocean referential) T-pt [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtau !: sea surface j-stress (ocean referential) T-pt [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utauU , utau_b !: sea surface i-stress (ocean referential) U-pt [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtauV , vtau_b !: sea surface j-stress (ocean referential) V-pt [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau_icb, vtau_icb !: sea surface (i,j)-stress used by icebergs [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: taum !: module of sea surface stress (at T-point) [N/m2]
!! wndm is used compute surface gases exchanges in ice-free ocean or leads
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wndm !: wind speed module at T-point (=|U10m-Uoce|) [m/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhoa !: air density at "rn_zu" m above the sea [kg/m3]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qsr !: sea heat flux: solar [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qns , qns_b !: sea heat flux: non solar [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qsr_tot !: total solar heat flux (over sea and ice) [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qns_tot !: total non solar heat flux (over sea and ice) [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: emp , emp_b !: freshwater budget: volume flux [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx , sfx_b !: salt flux [PSS.kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: emp_tot !: total E-P over ocean and ice [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fmmflx !: freshwater budget: freezing/melting [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rnf , rnf_b !: river runoff [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwficb , fwficb_b !: iceberg melting [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wndm !: wind speed module at T-point (=|U10m-Uoce|) [m/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhoa !: air density at "rn_zu" m above the sea [kg/m3]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qsr !: sea heat flux: solar [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qns , qns_b !: sea heat flux: non solar [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qsr_tot !: total solar heat flux (over sea and ice) [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qns_tot !: total non solar heat flux (over sea and ice) [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: emp , emp_b !: freshwater budget: volume flux [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx , sfx_b !: salt flux [PSS.kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: emp_tot !: total E-P over ocean and ice [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fmmflx !: freshwater budget: freezing/melting [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rnf , rnf_b !: river runoff [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwficb , fwficb_b !: iceberg melting [Kg/m2/s]
!!
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sbc_tsc, sbc_tsc_b !: sbc content trend [K.m/s] jpi,jpj,jpts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qsr_hc , qsr_hc_b !: heat content trend due to qsr flux [K.m/s] jpi,jpj,jpk
!!
!!
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tprecip !: total precipitation [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sprecip !: solid precipitation [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fr_i !: ice fraction = 1 - lead fraction (between 0 to 1)
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(:,:) :: tprecip !: total precipitation [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sprecip !: solid precipitation [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fr_i !: ice fraction = 1 - lead fraction (between 0 to 1)
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) [-]
!!---------------------------------------------------------------------
!! ABL Vertical Domain size
......@@ -177,8 +179,8 @@ CONTAINS
!!---------------------------------------------------------------------
ierr(:) = 0
!
ALLOCATE( utau(jpi,jpj) , utau_b(jpi,jpj) , taum(jpi,jpj) , &
& vtau(jpi,jpj) , vtau_b(jpi,jpj) , wndm(jpi,jpj) , rhoa(jpi,jpj) , STAT=ierr(1) )
ALLOCATE( utau(jpi,jpj) , utau_b(jpi,jpj) , utauU(jpi,jpj) , taum(jpi,jpj) , &
& vtau(jpi,jpj) , vtau_b(jpi,jpj) , vtauV(jpi,jpj) , wndm(jpi,jpj) , rhoa(jpi,jpj) , STAT=ierr(1) )
!
ALLOCATE( qns_tot(jpi,jpj) , qns (jpi,jpj) , qns_b(jpi,jpj), &
& qsr_tot(jpi,jpj) , qsr (jpi,jpj) , &
......@@ -205,9 +207,10 @@ CONTAINS
END FUNCTION sbc_oce_alloc
!!clem => this subroutine is never used in nemo
SUBROUTINE sbc_tau2wnd
!!---------------------------------------------------------------------
!! *** ROUTINE sbc_tau2wnd ***
!! *** ROUTINE ***
!!
!! ** Purpose : Estimation of wind speed as a function of wind stress
!!
......@@ -217,17 +220,14 @@ CONTAINS
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3
REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient
REAL(wp) :: ztx, zty, ztau, zcoef ! temporary variables
REAL(wp) :: ztau, zcoef ! temporary variables
INTEGER :: ji, jj ! dummy indices
!!---------------------------------------------------------------------
zcoef = 0.5 / ( zrhoa * zcdrag )
DO_2D( 0, 0, 0, 0 )
ztx = utau(ji-1,jj ) + utau(ji,jj)
zty = vtau(ji ,jj-1) + vtau(ji,jj)
ztau = SQRT( ztx * ztx + zty * zty )
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
ztau = SQRT( utau(ji,jj)*utau(ji,jj) + vtau(ji,jj)*vtau(ji,jj) )
wndm(ji,jj) = SQRT ( ztau * zcoef ) * tmask(ji,jj,1)
END_2D
CALL lbc_lnk( 'sbc_oce', wndm(:,:) , 'T', 1.0_wp )
!
END SUBROUTINE sbc_tau2wnd
......
src/OCE/SBC/sbcblk.F90
View file @ 78481193
......@@ -493,7 +493,7 @@ CONTAINS
!! the stress is assumed to be in the (i,j) mesh referential
!!
!! ** Action : defined at each time-step at the air-sea interface
!! - utau, vtau i- and j-component of the wind stress
!! - utau, vtau i- and j-component of the wind stress at T-point
!! - taum wind stress module at T-point
!! - wndm wind speed module at T-point over free ocean or leads in presence of sea-ice
!! - qns, qsr non-solar and solar heat fluxes
......@@ -611,10 +611,10 @@ CONTAINS
END SUBROUTINE sbc_blk
SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, pqair, & ! inp
& pslp , pst , pu , pv, & ! inp
& puatm, pvatm, pdqsr , pdqlw , & ! inp
& ptsk , pssq , pcd_du, psen, plat, pevp ) ! out
SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, pqair, & ! <<= in
& pslp , pst , pu , pv, & ! <<= in
& puatm, pvatm, pdqsr , pdqlw , & ! <<= in
& ptsk , pssq , pcd_du, psen, plat, pevp ) ! =>> out
!!---------------------------------------------------------------------
!! *** ROUTINE blk_oce_1 ***
!!
......@@ -657,7 +657,6 @@ CONTAINS
#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
REAL(wp), DIMENSION(jpi,jpj) :: zch_oce ! sensible heat transfert coefficient over ocean
......@@ -695,22 +694,20 @@ CONTAINS
#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) )
wndm(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
END_2D
#endif
! ----------------------------------------------------------------------------- !
! I Solar FLUX !
! ----------------------------------------------------------------------------- !
! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave
zztmp = 1. - albo
! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle
IF( ln_dm2dc ) THEN
qsr(:,:) = zztmp * sbc_dcy( pdqsr(:,:) ) * tmask(:,:,1)
qsr(:,:) = ( 1._wp - albo ) * sbc_dcy( pdqsr(:,:) ) * tmask(:,:,1)
ELSE
qsr(:,:) = zztmp * pdqsr(:,:) * tmask(:,:,1)
qsr(:,:) = ( 1._wp - albo ) * pdqsr(:,:) * tmask(:,:,1)
ENDIF
! ----------------------------------------------------------------------------- !
! II Turbulent FLUXES !
! ----------------------------------------------------------------------------- !
......@@ -718,69 +715,62 @@ CONTAINS
! specific humidity at SST
pssq(:,:) = rdct_qsat_salt * q_sat( ptsk(:,:), pslp(:,:) )
! Backup "bulk SST" and associated spec. hum.
IF( ln_skin_cs .OR. ln_skin_wl ) THEN
!! Backup "bulk SST" and associated spec. hum.
zztmp1(:,:) = zsspt(:,:)
zztmp2(:,:) = pssq(:,:)
zztmp2(:,:) = pssq (:,:)
ENDIF
!! Time to call the user-selected bulk parameterization for
!! == transfer coefficients ==! Cd, Ch, Ce at T-point, and more...
SELECT CASE( nblk )
! transfer coefficients (Cd, Ch, Ce at T-point, and more)
SELECT CASE( nblk ) ! user-selected bulk parameterization
!
CASE( np_NCAR )
CALL turb_ncar ( rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , &
& nb_iter=nn_iter_algo )
!
CASE( np_COARE_3p0 )
CALL turb_coare3p0( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
!
CASE( np_COARE_3p6 )
CALL turb_coare3p6( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
!
CASE( np_ECMWF )
CALL turb_ecmwf ( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
!
CASE( np_ANDREAS )
CALL turb_andreas ( rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo )
!
CASE DEFAULT
CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk parameterizaton selected' )
!
END SELECT
IF( iom_use('Cd_oce') ) CALL iom_put("Cd_oce", zcd_oce * tmask(:,:,1))
IF( iom_use('Ce_oce') ) CALL iom_put("Ce_oce", zce_oce * tmask(:,:,1))
IF( iom_use('Ch_oce') ) CALL iom_put("Ch_oce", zch_oce * tmask(:,:,1))
! outputs
IF( iom_use('Cd_oce') ) CALL iom_put( "Cd_oce", zcd_oce * tmask(:,:,1) )
IF( iom_use('Ce_oce') ) CALL iom_put( "Ce_oce", zce_oce * tmask(:,:,1) )
IF( iom_use('Ch_oce') ) CALL iom_put( "Ch_oce", zch_oce * tmask(:,:,1) )
!! LB: mainly here for debugging purpose:
IF( iom_use('theta_zt') ) CALL iom_put("theta_zt", (ptair-rt0) * tmask(:,:,1)) ! potential temperature at z=zt
IF( iom_use('q_zt') ) CALL iom_put("q_zt", pqair * tmask(:,:,1)) ! specific humidity "
IF( iom_use('theta_zu') ) CALL iom_put("theta_zu", (theta_zu -rt0) * tmask(:,:,1)) ! potential temperature at z=zu
IF( iom_use('q_zu') ) CALL iom_put("q_zu", q_zu * tmask(:,:,1)) ! specific humidity "
IF( iom_use('ssq') ) CALL iom_put("ssq", pssq * tmask(:,:,1)) ! saturation specific humidity at z=0
IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu * tmask(:,:,1)) ! bulk wind speed at z=zu
IF( iom_use('theta_zt') ) CALL iom_put( "theta_zt", (ptair-rt0) * tmask(:,:,1) ) ! potential temperature at z=zt
IF( iom_use('q_zt') ) CALL iom_put( "q_zt", pqair * tmask(:,:,1) ) ! specific humidity "
IF( iom_use('theta_zu') ) CALL iom_put( "theta_zu", (theta_zu -rt0) * tmask(:,:,1) ) ! potential temperature at z=zu
IF( iom_use('q_zu') ) CALL iom_put( "q_zu", q_zu * tmask(:,:,1) ) ! specific humidity "
IF( iom_use('ssq') ) CALL iom_put( "ssq", pssq * tmask(:,:,1) ) ! saturation specific humidity at z=0
IF( iom_use('wspd_blk') ) CALL iom_put( "wspd_blk", zU_zu * tmask(:,:,1) ) ! bulk wind speed at z=zu
! In the presence of sea-ice we do not use the cool-skin/warm-layer update of zsspt, pssq & ptsk from turb_*()
IF( ln_skin_cs .OR. ln_skin_wl ) THEN
!! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zsspt, pssq & ptsk:
WHERE ( fr_i(:,:) > 0.001_wp )
! sea-ice present, we forget about the update, using what we backed up before call to turb_*()
zsspt(:,:) = zztmp1(:,:)
pssq(:,:) = zztmp2(:,:)
pssq (:,:) = zztmp2(:,:)
END WHERE
! apply potential temperature increment to abolute SST
ptsk(:,:) = ptsk(:,:) + ( zsspt(:,:) - zztmp1(:,:) )
......@@ -809,10 +799,10 @@ CONTAINS
END_2D
CALL BULK_FORMULA( rn_zu, zsspt(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), &
& zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), &
& wndm(:,:), zU_zu(:,:), pslp(:,:), rhoa(:,:), &
& taum(:,:), psen(:,:), plat(:,:), &
& pEvap=pevp(:,:), pfact_evap=rn_efac )
& zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), &
& wndm(:,:), zU_zu(:,:), pslp(:,:), rhoa(:,:), &
& taum(:,:), psen(:,:), plat(:,:), &
& pEvap=pevp(:,:), pfact_evap=rn_efac )
psen(:,:) = psen(:,:) * tmask(:,:,1)
plat(:,:) = plat(:,:) * tmask(:,:,1)
......@@ -821,57 +811,42 @@ CONTAINS
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
IF( wndm(ji,jj) > 0._wp ) THEN
zztmp = taum(ji,jj) / wndm(ji,jj)
zztmp = taum(ji,jj) / wndm(ji,jj)
#if defined key_cyclone
ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
utau(ji,jj) = zztmp * zwnd_i(ji,jj)
vtau(ji,jj) = zztmp * zwnd_j(ji,jj)
#else
ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
utau(ji,jj) = zztmp * pwndi(ji,jj)
vtau(ji,jj) = zztmp * pwndj(ji,jj)
#endif
ELSE
ztau_i(ji,jj) = 0._wp
ztau_j(ji,jj) = 0._wp
utau(ji,jj) = 0._wp
vtau(ji,jj) = 0._wp
ENDIF
END_2D
IF( ln_crt_fbk ) THEN ! aply eq. 10 and 11 of Renault et al. 2020 (doi: 10.1029/2019MS001715)
zstmax = MIN( rn_stau_a * 3._wp + rn_stau_b, 0._wp ) ! set the max value of Stau corresponding to a wind of 3 m/s (<0)
DO_2D( 0, 1, 0, 1 ) ! end at jpj and jpi, as ztau_j(ji,jj+1) ztau_i(ji+1,jj) used in the next loop
zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax ) ! stau (<0) must be smaller than zstmax
ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj ) + pu(ji,jj) ) - puatm(ji,jj) )
ztau_j(ji,jj) = ztau_j(ji,jj) + zstau * ( 0.5_wp * ( pv(ji ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) )
taum(ji,jj) = SQRT( ztau_i(ji,jj) * ztau_i(ji,jj) + ztau_j(ji,jj) * ztau_j(ji,jj) )
DO_2D( 0, 0, 0, 0 )
zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax ) * tmask(ji,jj,1) ! stau (<0) must be smaller than zstmax
utau(ji,jj) = utau(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj ) + pu(ji,jj) ) - puatm(ji,jj) )
vtau(ji,jj) = vtau(ji,jj) + zstau * ( 0.5_wp * ( pv(ji ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) )
taum(ji,jj) = SQRT( utau(ji,jj) * utau(ji,jj) + vtau(ji,jj) * vtau(ji,jj) )
END_2D
CALL lbc_lnk( 'sbcblk', utau, 'T', -1._wp, vtau, 'T', -1._wp, taum, 'T', 1._wp )
ENDIF
! ... utau, vtau at U- and V_points, resp.
! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
! Note that coastal wind stress is not used in the code... so this extra care has no effect
DO_2D( 0, 0, 0, 0 ) ! start loop at 2, in case ln_crt_fbk = T
utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj ) ) &
& * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( ztau_j(ji,jj) + ztau_j(ji ,jj+1) ) &
& * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
END_2D
IF( ln_crt_fbk ) THEN
CALL lbc_lnk( 'sbcblk', utau, 'U', -1._wp, vtau, 'V', -1._wp, taum, 'T', 1._wp )
ELSE
CALL lbc_lnk( 'sbcblk', utau, 'U', -1._wp, vtau, 'V', -1._wp )
ENDIF
! Saving open-ocean wind-stress (module and components) on T-points:
CALL iom_put( "taum_oce", taum(:,:)*tmask(:,:,1) ) ! output wind stress module
!#LB: These 2 lines below mostly here for 'STATION_ASF' test-case, otherwize "utau" (U-grid) and vtau" (V-grid) does the job in: [DYN/dynatf.F90])
CALL iom_put( "utau_oce", ztau_i(:,:)*tmask(:,:,1) ) ! utau at T-points!
CALL iom_put( "vtau_oce", ztau_j(:,:)*tmask(:,:,1) ) ! vtau at T-points!
! Saving open-ocean wind-stress (module and components)
CALL iom_put( "taum_oce", taum(:,:) ) ! wind stress module
! ! LB: These 2 lines below mostly here for 'STATION_ASF' test-case
CALL iom_put( "utau_oce", utau(:,:) ) ! utau
CALL iom_put( "vtau_oce", vtau(:,:) ) ! vtau
IF(sn_cfctl%l_prtctl) THEN
CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ', mask1=tmask )
CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ', mask1=tmask )
CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, &
& tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask )
CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=tmask, &
& tab2d_2=vtau , clinfo2=' vtau : ', mask2=tmask )
CALL prt_ctl( tab2d_1=zcd_oce, clinfo1=' blk_oce_1: Cd : ', mask1=tmask )
ENDIF
!
......@@ -896,8 +871,8 @@ CONTAINS
!! at the ocean surface at each time step knowing Cd, Ch, Ce and
!! atmospheric variables (from ABL or external data)
!!
!! ** Outputs : - utau : i-component of the stress at U-point (N/m2)
!! - vtau : j-component of the stress at V-point (N/m2)
!! ** Outputs : - utau : i-component of the stress at T-point (N/m2)
!! - vtau : j-component of the stress at T-point (N/m2)
!! - taum : Wind stress module at T-point (N/m2)
!! - wndm : Wind speed module at T-point (m/s)
!! - qsr : Solar heat flux over the ocean (W/m2)
......@@ -1015,11 +990,16 @@ CONTAINS
REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pcd_dui ! if ln_abl
!
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zootm_su ! sea-ice surface mean temperature
REAL(wp) :: zztmp1, zztmp2 ! temporary scalars
REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zsipt ! temporary array
REAL(wp) :: zztmp ! temporary scalars
REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zsipt ! temporary array
REAL(wp), DIMENSION(jpi,jpj) :: zmsk00 ! O% concentration ice mask
!!---------------------------------------------------------------------
!
! treshold for outputs
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , fr_i(ji,jj) - 1.e-6_wp ) ) ! 1 if ice, 0 if no ice
END_2D
! ------------------------------------------------------------ !
! Wind module relative to the moving ice ( U10m - U_ice ) !
! ------------------------------------------------------------ !
......@@ -1032,9 +1012,9 @@ CONTAINS
zsipt(:,:) = theta_exner( ptsui(:,:), pslp(:,:) )
! sea-ice <-> atmosphere bulk transfer coefficients
SELECT CASE( nblk_ice )
CASE( np_ice_cst )
SELECT CASE( nblk_ice ) ! user-selected bulk parameterization
!
CASE( np_ice_cst )
! Constant bulk transfer coefficients over sea-ice:
Cd_ice(:,:) = rn_Cd_i
Ch_ice(:,:) = rn_Ch_i
......@@ -1042,62 +1022,45 @@ CONTAINS
! no height adjustment, keeping zt values:
theta_zu_i(:,:) = ptair(:,:)
q_zu_i(:,:) = pqair(:,:)
CASE( np_ice_an05 ) ! calculate new drag from Lupkes(2015) equations
!
CASE( np_ice_an05 ) ! from Andreas(2005) equations
ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
CALL turb_ice_an05( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice, &
& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
!!
CASE( np_ice_lu12 )
!
CASE( np_ice_lu12 ) ! from Lupkes(2012) equations
ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
CALL turb_ice_lu12( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice, fr_i, &
& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
!!
CASE( np_ice_lg15 ) ! calculate new drag from Lupkes(2015) equations
!
CASE( np_ice_lg15 ) ! from Lupkes and Gryanik (2015) equations
ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
CALL turb_ice_lg15( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice, fr_i, &
& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
!!
!
END SELECT
IF( iom_use('Cd_ice').OR.iom_use('Ce_ice').OR.iom_use('Ch_ice').OR.iom_use('taum_ice').OR.iom_use('utau_ice').OR.iom_use('vtau_ice') ) &
& ztmp(:,:) = ( 1._wp - MAX(0._wp, SIGN( 1._wp, 1.E-6_wp - fr_i )) )*tmask(:,:,1) ! mask for presence of ice !
IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice*ztmp)
IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice*ztmp)
IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice*ztmp)
IF( ln_blk ) THEN
! ---------------------------------------------------- !
! Wind stress relative to nonmoving ice ( U10m ) !
! ---------------------------------------------------- !
! supress moving ice in wind stress computation as we don't know how to do it properly...
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)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zztmp = rhoa(ji,jj) * Cd_ice(ji,jj) * wndm_ice(ji,jj)
putaui(ji,jj) = zztmp * pwndi(ji,jj)
pvtaui(ji,jj) = zztmp * pwndj(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 !!!
IF(iom_use('taum_ice')) CALL iom_put('taum_ice', SQRT( putaui*putaui + pvtaui*pvtaui )*ztmp )
!#LB: These 2 lines below mostly here for 'STATION_ASF' test-case, otherwize "utau_oi" (U-grid) and vtau_oi" (V-grid) does the job in: [ICE/icedyn_rhg_evp.F90])
IF(iom_use('utau_ice')) CALL iom_put("utau_ice", putaui*ztmp) ! utau at T-points!
IF(iom_use('vtau_ice')) CALL iom_put("vtau_ice", pvtaui*ztmp) ! vtau at T-points!
!
DO_2D( 0, 0, 0, 0 ) ! U & V-points (same as ocean).
!#LB: QUESTION?? so SI3 expects wind stress vector to be provided at U & V points? Not at T-points ?
! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology
zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) )
zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) )
putaui(ji,jj) = zztmp1 * ( putaui(ji,jj) + putaui(ji+1,jj ) )
pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji ,jj+1) )
END_2D
CALL lbc_lnk( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp )
! outputs
! LB: not weighted by the ice concentration
IF( iom_use('taum_ice') ) CALL iom_put( 'taum_ice', SQRT( putaui*putaui + pvtaui*pvtaui ) * zmsk00 )
! LB: These 2 lines below mostly here for 'STATION_ASF' test-case
IF( iom_use('utau_ice') ) CALL iom_put( "utau_ice", putaui * zmsk00 )
IF( iom_use('vtau_ice') ) CALL iom_put( "vtau_ice", pvtaui * zmsk00 )
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ', mask1=umask &
& , tab2d_2=pvtaui , clinfo2=' pvtaui : ', mask2=vmask )
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ', mask1=tmask &
& , tab2d_2=pvtaui , clinfo2=' pvtaui : ', mask2=tmask )
ELSE ! ln_abl
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
......@@ -1105,10 +1068,15 @@ CONTAINS
pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj)
pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj)
END_2D
pssqi(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ; ! more accurate way to obtain ssq !
pssqi(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! more accurate way to obtain ssq
ENDIF ! ln_blk / ln_abl
!
! outputs
IF( iom_use('Cd_ice') ) CALL iom_put( "Cd_ice", Cd_ice * zmsk00 )
IF( iom_use('Ce_ice') ) CALL iom_put( "Ce_ice", Ce_ice * zmsk00 )
IF( iom_use('Ch_ice') ) CALL iom_put( "Ch_ice", Ch_ice * zmsk00 )
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ', mask1=tmask )
!
END SUBROUTINE blk_ice_1
......
src/OCE/SBC/sbccpl.F90
View file @ 78481193
......@@ -68,15 +68,15 @@ MODULE sbccpl
INTEGER, PARAMETER :: jpr_otx1 = 1 ! 3 atmosphere-ocean stress components on grid 1
INTEGER, PARAMETER :: jpr_oty1 = 2 !
INTEGER, PARAMETER :: jpr_otz1 = 3 !
INTEGER, PARAMETER :: jpr_otx2 = 4 ! 3 atmosphere-ocean stress components on grid 2
INTEGER, PARAMETER :: jpr_oty2 = 5 !
INTEGER, PARAMETER :: jpr_otz2 = 6 !
!!$ INTEGER, PARAMETER :: jpr_otx2 = 4 ! 3 atmosphere-ocean stress components on grid 2
!!$ INTEGER, PARAMETER :: jpr_oty2 = 5 !
!!$ INTEGER, PARAMETER :: jpr_otz2 = 6 !
INTEGER, PARAMETER :: jpr_itx1 = 7 ! 3 atmosphere-ice stress components on grid 1
INTEGER, PARAMETER :: jpr_ity1 = 8 !
INTEGER, PARAMETER :: jpr_itz1 = 9 !
INTEGER, PARAMETER :: jpr_itx2 = 10 ! 3 atmosphere-ice stress components on grid 2
INTEGER, PARAMETER :: jpr_ity2 = 11 !
INTEGER, PARAMETER :: jpr_itz2 = 12 !
!!$ INTEGER, PARAMETER :: jpr_itx2 = 10 ! 3 atmosphere-ice stress components on grid 2
!!$ INTEGER, PARAMETER :: jpr_ity2 = 11 !
!!$ INTEGER, PARAMETER :: jpr_itz2 = 12 !
INTEGER, PARAMETER :: jpr_qsroce = 13 ! Qsr above the ocean
INTEGER, PARAMETER :: jpr_qsrice = 14 ! Qsr above the ice
INTEGER, PARAMETER :: jpr_qsrmix = 15
......@@ -128,9 +128,9 @@ MODULE sbccpl
INTEGER, PARAMETER :: jpr_isf = 60
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 :: jpr_qtrice = 63 ! Transmitted solar thru sea-ice
INTEGER, PARAMETER :: jprcv = 62 ! total number of fields received
INTEGER, PARAMETER :: jprcv = 63 ! total number of fields received
INTEGER, PARAMETER :: jps_fice = 1 ! ice fraction sent to the atmosphere
INTEGER, PARAMETER :: jps_toce = 2 ! ocean temperature
......@@ -194,7 +194,7 @@ MODULE sbccpl
& 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, &
& sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf, sn_rcv_ts_ice
& sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf, sn_rcv_ts_ice, sn_rcv_qtrice
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
......@@ -202,7 +202,6 @@ MODULE sbccpl
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
! ! 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)
......@@ -281,7 +280,7 @@ CONTAINS
& 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_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_iceflx, sn_rcv_co2 , sn_rcv_icb , sn_rcv_isf , sn_rcv_ts_ice, sn_rcv_qtrice, &
& sn_rcv_mslp
!!---------------------------------------------------------------------
......@@ -313,7 +312,6 @@ CONTAINS
WRITE(numout,*)' surface stress = ', TRIM(sn_rcv_tau%cldes ), ' (', TRIM(sn_rcv_tau%clcat ), ')'
WRITE(numout,*)' - referential = ', sn_rcv_tau%clvref
WRITE(numout,*)' - orientation = ', sn_rcv_tau%clvor
WRITE(numout,*)' - mesh = ', sn_rcv_tau%clvgrd
WRITE(numout,*)' non-solar heat flux sensitivity = ', TRIM(sn_rcv_dqnsdt%cldes), ' (', TRIM(sn_rcv_dqnsdt%clcat), ')'
WRITE(numout,*)' solar heat flux = ', TRIM(sn_rcv_qsr%cldes ), ' (', TRIM(sn_rcv_qsr%clcat ), ')'
WRITE(numout,*)' non-solar heat flux = ', TRIM(sn_rcv_qns%cldes ), ' (', TRIM(sn_rcv_qns%clcat ), ')'
......@@ -323,7 +321,7 @@ CONTAINS
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 thru sea-ice = ', TRIM(sn_rcv_qtrice%cldes), ' (', TRIM(sn_rcv_qtrice%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:'
......@@ -358,7 +356,7 @@ CONTAINS
WRITE(numout,*)' - mesh = ', sn_snd_crtw%clvgrd
ENDIF
IF( lwp .AND. ln_wave) THEN ! control print
WRITE(numout,*)' surface waves:'
WRITE(numout,*)' surface waves:'
WRITE(numout,*)' Significant wave heigth = ', TRIM(sn_rcv_hsig%cldes ), ' (', TRIM(sn_rcv_hsig%clcat ), ')'
WRITE(numout,*)' Wave to oce energy flux = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')'
WRITE(numout,*)' Surface Stokes drift grid u = ', TRIM(sn_rcv_sdrfx%cldes ), ' (', TRIM(sn_rcv_sdrfx%clcat ), ')'
......@@ -368,8 +366,8 @@ CONTAINS
WRITE(numout,*)' Stress frac adsorbed by waves = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')'
WRITE(numout,*)' Neutral surf drag coefficient = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')'
WRITE(numout,*)' Charnock coefficient = ', TRIM(sn_rcv_charn%cldes ), ' (', TRIM(sn_rcv_charn%clcat ), ')'
WRITE(numout,*)' Transport associated to Stokes drift grid u = ', TRIM(sn_rcv_tusd%cldes ), ' (', TRIM(sn_rcv_tusd%clcat ), ')'
WRITE(numout,*)' Transport associated to Stokes drift grid v = ', TRIM(sn_rcv_tvsd%cldes ), ' (', TRIM(sn_rcv_tvsd%clcat ), ')'
WRITE(numout,*)' Transport associated to Stokes drift u = ', TRIM(sn_rcv_tusd%cldes ), ' (', TRIM(sn_rcv_tusd%clcat ), ')'
WRITE(numout,*)' Transport associated to Stokes drift v = ', TRIM(sn_rcv_tvsd%cldes ), ' (', TRIM(sn_rcv_tvsd%clcat ), ')'
WRITE(numout,*)' Bernouilli pressure head = ', TRIM(sn_rcv_bhd%cldes ), ' (', TRIM(sn_rcv_bhd%clcat ), ')'
WRITE(numout,*)'Wave to ocean momentum flux and Net wave-supported stress = ', TRIM(sn_rcv_taw%cldes ), ' (', TRIM(sn_rcv_taw%clcat ), ')'
WRITE(numout,*)' Surface current to waves = ', TRIM(sn_snd_crtw%cldes ), ' (', TRIM(sn_snd_crtw%clcat ), ')'
......@@ -399,87 +397,26 @@ CONTAINS
srcv(jpr_otx1)%clname = 'O_OTaux1' ! 1st ocean component on grid ONE (T or U)
srcv(jpr_oty1)%clname = 'O_OTauy1' ! 2nd - - - -
srcv(jpr_otz1)%clname = 'O_OTauz1' ! 3rd - - - -
srcv(jpr_otx2)%clname = 'O_OTaux2' ! 1st ocean component on grid TWO (V)
srcv(jpr_oty2)%clname = 'O_OTauy2' ! 2nd - - - -
srcv(jpr_otz2)%clname = 'O_OTauz2' ! 3rd - - - -
!
srcv(jpr_itx1)%clname = 'O_ITaux1' ! 1st ice component on grid ONE (T, F, I or U)
srcv(jpr_ity1)%clname = 'O_ITauy1' ! 2nd - - - -
srcv(jpr_itz1)%clname = 'O_ITauz1' ! 3rd - - - -
srcv(jpr_itx2)%clname = 'O_ITaux2' ! 1st ice component on grid TWO (V)
srcv(jpr_ity2)%clname = 'O_ITauy2' ! 2nd - - - -
srcv(jpr_itz2)%clname = 'O_ITauz2' ! 3rd - - - -
!
! Vectors: change of sign at north fold ONLY if on the local grid
IF( TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM( sn_rcv_tau%cldes ) == 'oce and ice' &
.OR. TRIM( sn_rcv_tau%cldes ) == 'mixed oce-ice' ) THEN ! avoid working with the atmospheric fields if they are not coupled
!
IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' ) srcv(jpr_otx1:jpr_itz2)%nsgn = -1.
! ! Set grid and action
SELECT CASE( TRIM( sn_rcv_tau%clvgrd ) ) ! 'T', 'U,V', 'U,V,I', 'U,V,F', 'T,I', 'T,F', or 'T,U,V'
CASE( 'T' )
srcv(jpr_otx1:jpr_itz2)%clgrid = 'T' ! oce and ice components given at T-point
srcv(jpr_otx1:jpr_otz1)%laction = .TRUE. ! receive oce components on grid 1
srcv(jpr_itx1:jpr_itz1)%laction = .TRUE. ! receive ice components on grid 1
CASE( 'U,V' )
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 = 'U' ! ice components given at U-point
srcv(jpr_itx2:jpr_itz2)%clgrid = 'V' ! and V-point
srcv(jpr_otx1:jpr_itz2)%laction = .TRUE. ! receive oce and ice components on both grid 1 & 2
CASE( 'U,V,T' )
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 = 'T' ! ice components given at T-point
srcv(jpr_otx1:jpr_otz2)%laction = .TRUE. ! receive oce components on grid 1 & 2
srcv(jpr_itx1:jpr_itz1)%laction = .TRUE. ! receive ice components on grid 1 only
CASE( 'U,V,I' )
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 = 'I' ! ice components given at I-point
srcv(jpr_otx1:jpr_otz2)%laction = .TRUE. ! receive oce components on grid 1 & 2
srcv(jpr_itx1:jpr_itz1)%laction = .TRUE. ! receive ice components on grid 1 only
CASE( 'U,V,F' )
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_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
srcv(jpr_itx1:jpr_itz1)%clgrid = 'I' ! ice components given at I-point
srcv(jpr_otx1:jpr_otz1)%laction = .TRUE. ! receive oce components on grid 1
srcv(jpr_itx1:jpr_itz1)%laction = .TRUE. ! receive ice components on grid 1
CASE( 'T,F' )
srcv(jpr_otx1:jpr_itz2)%clgrid = 'T' ! oce and ice components given at T-point
srcv(jpr_itx1:jpr_itz1)%clgrid = 'F' ! ice components given at F-point
srcv(jpr_otx1:jpr_otz1)%laction = .TRUE. ! receive oce components on grid 1
srcv(jpr_itx1:jpr_itz1)%laction = .TRUE. ! receive ice components on grid 1
CASE( 'T,U,V' )
srcv(jpr_otx1:jpr_otz1)%clgrid = 'T' ! oce components given at T-point
srcv(jpr_itx1:jpr_itz1)%clgrid = 'U' ! ice components given at U-point
srcv(jpr_itx2:jpr_itz2)%clgrid = 'V' ! and V-point
srcv(jpr_otx1:jpr_otz1)%laction = .TRUE. ! receive oce components on grid 1 only
srcv(jpr_itx1:jpr_itz2)%laction = .TRUE. ! receive ice components on grid 1 & 2
CASE default
CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_tau%clvgrd' )
END SELECT
IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' ) srcv(jpr_otx1:jpr_itz1)%nsgn = -1.
! ! Set grid and action
srcv(jpr_otx1:jpr_otz1)%laction = .TRUE. ! receive oce components on grid 1
srcv(jpr_itx1:jpr_itz1)%laction = .TRUE. ! receive ice components on grid 1
!
IF( TRIM( sn_rcv_tau%clvref ) == 'spherical' ) & ! spherical: 3rd component not received
& srcv( (/jpr_otz1, jpr_otz2, jpr_itz1, jpr_itz2/) )%laction = .FALSE.
!
IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' ) THEN ! already on local grid -> no need of the second grid
srcv(jpr_otx2:jpr_otz2)%laction = .FALSE.
srcv(jpr_itx2:jpr_itz2)%laction = .FALSE.
srcv(jpr_oty1)%clgrid = srcv(jpr_oty2)%clgrid ! not needed but cleaner...
srcv(jpr_ity1)%clgrid = srcv(jpr_ity2)%clgrid ! not needed but cleaner...
ENDIF
& srcv( (/jpr_otz1, jpr_itz1/) )%laction = .FALSE.
!
IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN ! 'oce and ice' case ocean stress on ocean mesh used
srcv(jpr_itx1:jpr_itz2)%laction = .FALSE. ! ice components not received
srcv(jpr_itx1)%clgrid = 'U' ! ocean stress used after its transformation
srcv(jpr_ity1)%clgrid = 'V' ! i.e. it is always at U- & V-points for i- & j-comp. resp.
IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN ! 'oce and ice' case ocean stress on ocean mesh used
srcv(jpr_itx1:jpr_itz1)%laction = .FALSE. ! ice components not received
ENDIF
ENDIF
......@@ -612,18 +549,18 @@ CONTAINS
ENDIF
srcv(jpr_topm:jpr_botm)%laction = .TRUE.
ENDIF
!! ! ! --------------------------- !
!! ! ! transmitted solar thru ice !
!! ! ! --------------------------- !
!! srcv(jpr_qtrice)%clname = 'OQtr'
!! IF( TRIM(sn_rcv_qtrice%cldes) == 'coupled' ) THEN
!! IF ( TRIM( sn_rcv_qtrice%clcat ) == 'yes' ) THEN
!! srcv(jpr_qtrice)%nct = nn_cats_cpl
!! ELSE
!! CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qtrice%clcat should always be set to yes currently' )
!! ENDIF
!! srcv(jpr_qtrice)%laction = .TRUE.
!! ENDIF
! ! --------------------------- !
! ! transmitted solar thru ice !
! ! --------------------------- !
srcv(jpr_qtrice)%clname = 'OQtr'
IF( TRIM(sn_rcv_qtrice%cldes) == 'coupled' ) THEN
IF ( TRIM( sn_rcv_qtrice%clcat ) == 'yes' ) THEN
srcv(jpr_qtrice)%nct = nn_cats_cpl
ELSE
CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qtrice%clcat should always be set to yes currently' )
ENDIF
srcv(jpr_qtrice)%laction = .TRUE.
ENDIF
! ! ------------------------- !
! ! ice skin temperature !
! ! ------------------------- !
......@@ -725,11 +662,10 @@ CONTAINS
srcv(:)%clgrid = 'T' ! force default definition in case of opa <-> sas coupling
srcv(:)%nsgn = 1. ! force default definition in case of opa <-> sas coupling
srcv( (/jpr_qsroce, jpr_qnsoce, jpr_oemp, jpr_sflx, jpr_fice, jpr_otx1, jpr_oty1, jpr_taum/) )%laction = .TRUE.
srcv(jpr_otx1)%clgrid = 'U' ! oce components given at U-point
srcv(jpr_oty1)%clgrid = 'V' ! and V-point
srcv(jpr_otx1)%clgrid = 'T' ! oce components given at T-point
srcv(jpr_oty1)%clgrid = 'T'
! Vectors: change of sign at north fold ONLY if on the local grid
srcv( (/jpr_otx1,jpr_oty1/) )%nsgn = -1.
sn_rcv_tau%clvgrd = 'U,V'
sn_rcv_tau%clvor = 'local grid'
sn_rcv_tau%clvref = 'spherical'
sn_rcv_emp%cldes = 'oce only'
......@@ -1162,7 +1098,7 @@ CONTAINS
!! ** Method : receive all fields from the atmosphere and transform
!! them into ocean surface boundary condition fields
!!
!! ** Action : update utau, vtau ocean stress at U,V grid
!! ** Action : update utau, vtau ocean stress at T-point
!! taum wind stress module at T-point
!! wndm wind speed module at T-point over free ocean or leads in presence of sea-ice
!! qns non solar heat fluxes including emp heat content (ocean only case)
......@@ -1221,39 +1157,20 @@ CONTAINS
IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN ! 2 components on the sphere
! ! (cartesian to spherical -> 3 to 2 components)
!
CALL geo2oce( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), frcv(jpr_otz1)%z3(:,:,1), &
& srcv(jpr_otx1)%clgrid, ztx, zty )
CALL geo2oce( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), frcv(jpr_otz1)%z3(:,:,1), 'T', ztx, zty )
frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:) ! overwrite 1st comp. on the 1st grid
frcv(jpr_oty1)%z3(:,:,1) = zty(:,:) ! overwrite 2nd comp. on the 1st grid
!
IF( srcv(jpr_otx2)%laction ) THEN
CALL geo2oce( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), frcv(jpr_otz2)%z3(:,:,1), &
& srcv(jpr_otx2)%clgrid, ztx, zty )
frcv(jpr_otx2)%z3(:,:,1) = ztx(:,:) ! overwrite 1st comp. on the 2nd grid
frcv(jpr_oty2)%z3(:,:,1) = zty(:,:) ! overwrite 2nd comp. on the 2nd grid
ENDIF
!
ENDIF
!
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 )
ELSE
CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty )
ENDIF
CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), 'T', 'en->i', ztx )
CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), 'T', 'en->j', zty )
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
DO_2D( 0, 0, 0, 0 ) ! T ==> (U,V)
frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj ,1) + frcv(jpr_otx1)%z3(ji,jj,1) )
frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) )
END_2D
CALL lbc_lnk( 'sbccpl', frcv(jpr_otx1)%z3(:,:,1), 'U', -1.0_wp, frcv(jpr_oty1)%z3(:,:,1), 'V', -1.0_wp )
ENDIF
llnewtx = .TRUE.
ELSE
llnewtx = .FALSE.
......@@ -1272,12 +1189,9 @@ CONTAINS
IF( .NOT. srcv(jpr_taum)%laction ) THEN ! compute wind stress module from its components if not received
! => need to be done only when otx1 was changed
IF( llnewtx ) THEN
DO_2D( 0, 0, 0, 0 )
zzx = frcv(jpr_otx1)%z3(ji-1,jj ,1) + frcv(jpr_otx1)%z3(ji,jj,1)
zzy = frcv(jpr_oty1)%z3(ji ,jj-1,1) + frcv(jpr_oty1)%z3(ji,jj,1)
frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
END_2D
CALL lbc_lnk( 'sbccpl', frcv(jpr_taum)%z3(:,:,1), 'T', 1.0_wp )
zzx = frcv(jpr_otx1)%z3(ji,jj,1)
zzy = frcv(jpr_oty1)%z3(ji,jj,1)
frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
llnewtau = .TRUE.
ELSE
llnewtau = .FALSE.
......@@ -1594,7 +1508,7 @@ CONTAINS
!! ** Action : return ptau_i, ptau_j, the stress over the ice
!!----------------------------------------------------------------------
REAL(wp), INTENT(inout), DIMENSION(:,:) :: p_taui ! i- & j-components of atmos-ice stress [N/m2]
REAL(wp), INTENT(inout), DIMENSION(:,:) :: p_tauj ! at I-point (B-grid) or U & V-point (C-grid)
REAL(wp), INTENT(inout), DIMENSION(:,:) :: p_tauj ! at T-point
!!
INTEGER :: ji, jj ! dummy loop indices
INTEGER :: itx ! index of taux over ice
......@@ -1616,28 +1530,16 @@ CONTAINS
!
IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN ! 2 components on the sphere
! ! (cartesian to spherical -> 3 to 2 components)
CALL geo2oce( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), frcv(jpr_itz1)%z3(:,:,1), &
& srcv(jpr_itx1)%clgrid, ztx, zty )
CALL geo2oce( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), frcv(jpr_itz1)%z3(:,:,1), 'T', ztx, zty )
frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:) ! overwrite 1st comp. on the 1st grid
frcv(jpr_ity1)%z3(:,:,1) = zty(:,:) ! overwrite 2nd comp. on the 1st grid
!
IF( srcv(jpr_itx2)%laction ) THEN
CALL geo2oce( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), frcv(jpr_itz2)%z3(:,:,1), &
& srcv(jpr_itx2)%clgrid, ztx, zty )
frcv(jpr_itx2)%z3(:,:,1) = ztx(:,:) ! overwrite 1st comp. on the 2nd grid
frcv(jpr_ity2)%z3(:,:,1) = zty(:,:) ! overwrite 2nd comp. on the 2nd grid
ENDIF
!
ENDIF
!
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_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->i', ztx )
IF( srcv(jpr_itx2)%laction ) THEN
CALL rot_rep( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), srcv(jpr_itx2)%clgrid, 'en->j', zty )
ELSE
CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->j', zty )
ENDIF
CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), 'T', 'en->i', ztx )
CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), 'T', 'en->j', zty )
frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:) ! overwrite 1st component on the 1st grid
frcv(jpr_ity1)%z3(:,:,1) = zty(:,:) ! overwrite 2nd component on the 1st grid
ENDIF
......@@ -1651,29 +1553,8 @@ CONTAINS
! ! ======================= !
! ! put on ice grid !
! ! ======================= !
!
! j+1 j -----V---F
! ice stress on ice velocity point ! |
! (C-grid ==>(U,V)) j | T U
! | |
! j j-1 -I-------|
! (for I) | |
! i-1 i i
! i i+1 (for I)
SELECT CASE ( srcv(jpr_itx1)%clgrid )
CASE( 'U' )
p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1) ! (U,V) ==> (U,V)
p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
CASE( 'T' )
DO_2D( 0, 0, 0, 0 ) ! T ==> (U,V)
! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology
zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) )
zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) )
p_taui(ji,jj) = zztmp1 * ( frcv(jpr_itx1)%z3(ji+1,jj ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
p_tauj(ji,jj) = zztmp2 * ( frcv(jpr_ity1)%z3(ji ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
END_2D
CALL lbc_lnk( 'sbccpl', p_taui, 'U', -1., p_tauj, 'V', -1. )
END SELECT
p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)
p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
ENDIF
!
......@@ -1794,8 +1675,8 @@ CONTAINS
END SELECT
! --- evaporation over ice (kg/m2/s) --- !
IF (ln_scale_ice_flux) THEN ! typically met-office requirements
IF (sn_rcv_emp%clcat == 'yes') THEN
IF( ln_scale_ice_flux ) THEN ! typically met-office requirements
IF( sn_rcv_emp%clcat == 'yes' ) THEN
WHERE( a_i(:,:,:) > 1.e-10 ) ; zevap_ice(:,:,:) = frcv(jpr_ievp)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
ELSEWHERE ; zevap_ice(:,:,:) = 0._wp
END WHERE
......@@ -1812,7 +1693,7 @@ CONTAINS
ENDDO
ENDIF
ELSE
IF (sn_rcv_emp%clcat == 'yes') THEN
IF( sn_rcv_emp%clcat == 'yes' ) THEN
zevap_ice(:,:,1:jpl) = frcv(jpr_ievp)%z3(:,:,1:jpl)
WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:)
ELSEWHERE ; zevap_ice_total(:,:) = 0._wp
......@@ -1826,7 +1707,7 @@ CONTAINS
ENDIF
ENDIF
IF ( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN
IF( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN
! For conservative case zemp_ice has not been defined yet. Do it now.
zemp_ice(:,:) = zevap_ice_total(:,:) * picefr(:,:) - frcv(jpr_snow)%z3(:,:,1) * picefr(:,:)
ENDIF
......@@ -1926,13 +1807,13 @@ CONTAINS
& - zevap_ice_total(:,:) * picefr(:,:) ) * tmask(:,:,1) ) ! ice-free oce evap (cell average)
! note: runoff output is done in sbcrnf (which includes icebergs too) and iceshelf output is done in sbcisf
!! 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
!! 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
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(:,:,:)
......@@ -2213,28 +2094,29 @@ 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_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
!
ENDIF
......@@ -2403,10 +2285,10 @@ CONTAINS
CASE( 'weighted ice' ) ;
SELECT CASE( sn_snd_alb%clcat )
CASE( 'yes' )
ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
CASE( 'no' )
WHERE( fr_i (:,:) > 0. )
ztmp1(:,:) = SUM ( alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 )
ztmp1(:,:) = SUM ( alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 )
ELSEWHERE
ztmp1(:,:) = 0.
END WHERE
......@@ -2572,17 +2454,18 @@ CONTAINS
IF( ssnd(jps_ocx1)%laction ) THEN ! Surface current !
! ! ------------------------- !
!
! j+1 j -----V---F
! surface velocity always sent from T point ! |
! j | T U
! | |
! j j-1 -I-------|
! (for I) | |
! i-1 i i
! i i+1 (for I)
! j -----V---F
! surface velocity always sent from T point ! |
! j | T U
! | |
! j-1 -I-------|
! | |
! i-1 i i
!!clem: make a new variable at T-point to replace uu and vv => uuT and vvT for instance
IF( nn_components == jp_iam_oce ) THEN
zotx1(:,:) = uu(:,:,1,Kmm)
zoty1(:,:) = vv(:,:,1,Kmm)
!!clem : should be demi sum, no? Or uuT and vvT
ELSE
SELECT CASE( TRIM( sn_snd_crt%cldes ) )
CASE( 'oce only' ) ! C-grid ==> T
......@@ -2652,42 +2535,42 @@ CONTAINS
! ! Surface current to waves !
! ! ------------------------- !
IF( ssnd(jps_ocxw)%laction .OR. ssnd(jps_ocyw)%laction ) THEN
!
! j+1 j -----V---F
! surface velocity always sent from T point ! |
! j | T U
! | |
! j j-1 -I-------|
! (for I) | |
! i-1 i i
! i i+1 (for I)
SELECT CASE( TRIM( sn_snd_crtw%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
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)
zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * fr_i(ji,jj)
zity1(ji,jj) = 0.5 * ( v_ice(ji,jj ) + v_ice(ji ,jj-1 ) ) * fr_i(ji,jj)
END_2D
CALL lbc_lnk( 'sbccpl', zitx1, 'T', -1.0_wp, zity1, 'T', -1.0_wp )
CASE( 'mixed oce-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) &
& + 0.5 * ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * fr_i(ji,jj)
zoty1(ji,jj) = 0.5 * ( vv (ji,jj,1,Kmm) + vv (ji ,jj-1,1,Kmm) ) * zfr_l(ji,jj) &
& + 0.5 * ( v_ice(ji,jj ) + v_ice(ji ,jj-1 ) ) * fr_i(ji,jj)
END_2D
END SELECT
!
! j -----V---F
! surface velocity always sent from T point ! |
! j | T U
! | |
! j-1 -I-------|
! | |
! i-1 i i
!!clem: make a new variable at T-point to replace uu and vv => uuT and vvT for instance
SELECT CASE( TRIM( sn_snd_crtw%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
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)
zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * fr_i(ji,jj)
zity1(ji,jj) = 0.5 * ( v_ice(ji,jj ) + v_ice(ji ,jj-1 ) ) * fr_i(ji,jj)
END_2D
CALL lbc_lnk( 'sbccpl', zitx1, 'T', -1.0_wp, zity1, 'T', -1.0_wp )
CASE( 'mixed oce-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) &
& + 0.5 * ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * fr_i(ji,jj)
zoty1(ji,jj) = 0.5 * ( vv (ji,jj,1,Kmm) + vv (ji ,jj-1,1,Kmm) ) * zfr_l(ji,jj) &
& + 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_ocxw)%clgrid, -1.0_wp, zoty1, ssnd(jps_ocyw)%clgrid, -1.0_wp )
!
!
IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) THEN ! Rotation of the components
! ! Ocean component
! ! Ocean component
CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->e', ztmp1 ) ! 1st component
CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->n', ztmp2 ) ! 2nd component
zotx1(:,:) = ztmp1(:,:) ! overwrite the components
......@@ -2700,18 +2583,18 @@ CONTAINS
ENDIF
ENDIF
!
! ! spherical coordinates to cartesian -> 2 components to 3 components
! IF( TRIM( sn_snd_crtw%clvref ) == 'cartesian' ) THEN
! ztmp1(:,:) = zotx1(:,:) ! ocean currents
! ztmp2(:,:) = zoty1(:,:)
! CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
! !
! IF( ssnd(jps_ivx1)%laction ) THEN ! ice velocities
! ztmp1(:,:) = zitx1(:,:)
! ztmp1(:,:) = zity1(:,:)
! CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
! ENDIF
! ENDIF
! ! spherical coordinates to cartesian -> 2 components to 3 components
! IF( TRIM( sn_snd_crtw%clvref ) == 'cartesian' ) THEN
! ztmp1(:,:) = zotx1(:,:) ! ocean currents
! ztmp2(:,:) = zoty1(:,:)
! CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
! !
! IF( ssnd(jps_ivx1)%laction ) THEN ! ice velocities
! ztmp1(:,:) = zitx1(:,:)
! ztmp1(:,:) = zity1(:,:)
! CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
! ENDIF
! ENDIF
!
IF( ssnd(jps_ocxw)%laction ) CALL cpl_snd( jps_ocxw, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info ) ! ocean x current 1st grid
IF( ssnd(jps_ocyw)%laction ) CALL cpl_snd( jps_ocyw, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info ) ! ocean y current 1st grid
......