MODULE sbccpl !!====================================================================== !! *** MODULE sbccpl *** !! Ocean forcing: momentum, heat and freshwater coupled formulation !!===================================================================== !! History : 9.0 ! 06-07 (R. Redler, N. Keenlyside, W. Park) !! Original code split into flxmod & taumod !! 9.0 ! 06-07 (G. Madec) surface module !!---------------------------------------------------------------------- #if defined key_sbc_cpl !!---------------------------------------------------------------------- !! 'key_sbc_cpl' Coupled Ocean/Atmosphere formulation !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! namsbc_cpl : coupled formulation namlist !! sbc_cpl : coupled formulation for the ocean surface boundary condition !!---------------------------------------------------------------------- USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE phycst ! physical constants USE in_out_manager ! I/O manager USE lib_mpp ! distribued memory computing library USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE daymod ! calendar USE cpl_oasis3 ! OASIS3 coupling (to ECHAM5) USE cpl_oasis4 ! OASIS4 coupling (to ECHAM5) USE geo2ocean, ONLY : repere, repcmo USE ice_2, only : frld ! : leads fraction = 1-a/totalarea USE sbc_oce ! Surface boundary condition: ocean fields USE iom ! NetCDF library IMPLICIT NONE PRIVATE PUBLIC sbc_cpl ! routine called by step.F90 LOGICAL, PUBLIC :: lk_sbc_cpl = .TRUE. !: coupled formulation flag INTEGER , PARAMETER :: jpfld = 5 ! maximum number of files to read INTEGER , PARAMETER :: jp_taux = 1 ! index of wind stress (i-component) file INTEGER , PARAMETER :: jp_tauy = 2 ! index of wind stress (j-component) file INTEGER , PARAMETER :: jp_qtot = 3 ! index of total (non solar+solar) heat file INTEGER , PARAMETER :: jp_qsr = 4 ! index of solar heat file INTEGER , PARAMETER :: jp_emp = 5 ! index of evaporation-precipation file !!wonsun REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & taux, tauy & !: surface stress components in (i,j) referential USE sbc_ice, only : dqns_ice , & ! : derivative of non solar heat flux on sea ice qsr_ice , & ! : solar flux over ice qns_ice , & ! : total non solar heat flux (Longwave downward radiation) over ice tn_ice , & ! : ice surface temperature alb_ice , & ! : albedo of ice sprecip , & ! : solid (snow) precipitation over water (!) what about ice? tprecip , & ! : total precipitation ( or liquid precip minus evaporation in coupled mode) calving , & ! : calving rrunoff , & ! : monthly runoff (kg/m2/s) fr1_i0 , & ! : 1st part of the fraction of sol.rad. which penetrate inside the ice cover fr2_i0 ! : 2nd part of the fraction of sol.rad. which penetrate inside the ice cover USE ice_2, only : hicif , & ! : ice thickness frld , & ! : leads fraction = 1-a/totalarea hsnif , & ! : snow thickness u_ice , v_ice ! : ice velocity USE sbc_oce, only : sst_m ! : sea surface temperature REAL(wp), PUBLIC :: & !!! surface fluxes namelist (namflx) q0 = 0.e0, & ! net heat flux qsr0 = 0.e0, & ! solar heat flux emp0 = 0.e0, & ! net freshwater flux dqdt0 = -40., & ! coefficient for SST damping (W/m2/K) deds0 = 27.7 ! coefficient for SSS damping (mm/day) REAL(wp), DIMENSION(jpi,jpj) :: qsr_oce_recv , qsr_ice_recv REAL(wp), DIMENSION(jpi,jpj) :: qns_oce_recv, qns_ice_recv REAL(wp), DIMENSION(jpi,jpj) :: dqns_ice_recv REAL(wp), DIMENSION(jpi,jpj) :: tprecip_recv , precip_recv REAL(wp), DIMENSION(jpi,jpj) :: fr1_i0_recv , fr2_i0_recv REAL(wp), DIMENSION(jpi,jpj) :: rrunoff_recv , calving_recv #if defined key_cpl_ocevel REAL(wp), DIMENSION(jpi,jpj) :: un_weighted, vn_weighted REAL(wp), DIMENSION(jpi,jpj) :: un_send , vn_send #endif REAL(wp), DIMENSION(jpi,jpj) :: zrunriv ! river discharge into ocean REAL(wp), DIMENSION(jpi,jpj) :: zruncot ! continental discharge into ocean REAL(wp), DIMENSION(jpi,jpj) :: zpew ! P-E over water REAL(wp), DIMENSION(jpi,jpj) :: zpei ! P-E over ice REAL(wp), DIMENSION(jpi,jpj) :: zpsol ! surface downward snow fall REAL(wp), DIMENSION(jpi,jpj) :: zevice ! surface upward snow flux where sea ice !!wonsun !! * Substitutions # include "domzgr_substitute.h90" !!---------------------------------------------------------------------- !! OPA 9.0 , LOCEAN-IPSL (2006) !! $ Id: $ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE sbc_cpl( kt ) !!--------------------------------------------------------------------- !! *** ROUTINE sbc_cpl *** !! !! ** Purpose : provide at each time step the surface ocean fluxes !! (momentum, heat, freshwater and runoff) in coupled mode !! !! ** Method : - Recieve from a Atmospheric model via OASIS coupler : !! i-component of the stress taux (N/m2) !! j-component of the stress tauy (N/m2) !! net downward heat flux qtot (watt/m2) !! net downward radiative flux qsr (watt/m2) !! net upward freshwater (evapo - precip) emp (kg/m2/s) !! - send to the Atmospheric model via OASIS coupler : !! !! ** Action : update at each time-step the two components of the !! surface stress in both (i,j) and geographical ref. !! !! !! CAUTION : - never mask the surface stress fields !! !! ** Action : update at each time-step !! - taux & tauy : stress components in (i,j) referential !! - qns : non solar heat flux !! - qsr : solar heat flux !! - emp : evap - precip (volume flux) !! - emps : evap - precip (concentration/dillution) !! !! References : The OASIS User Guide, Version 3.0 and 4.0 !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: kt ! ocean time step !! INTEGER :: ji, jj ! dummy loop indices #if defined key_cpl_ocevel INTEGER :: ikchoix #endif INTEGER :: var_id, info INTEGER :: date !???? !!gm bug this is a real !!! REAL(wp) :: zfacflx, zfacwat, zfact REAL(wp), DIMENSION(jpi,jpj) :: ztaueuw, ztauevw ! eastward wind stress over water at U and V-points REAL(wp), DIMENSION(jpi,jpj) :: ztaunuw, ztaunvw ! northward wind stress over water at U and V-points REAL(wp), DIMENSION(jpi,jpj) :: ztaueui, ztauevi ! eastward wind stress over ice at U and V-points REAL(wp), DIMENSION(jpi,jpj) :: ztaunui, ztaunvi ! northward wind stress over ice at U and V-points REAL(wp), DIMENSION(jpi,jpj) :: ztaueu , ztauev ! eastward wind stress combined REAL(wp), DIMENSION(jpi,jpj) :: ztaunu , ztaunv ! northward wind stress combined !!--------------------------------------------------------------------- date = ( kt - nit000 ) * rdttra(1) ! date of exxhanges ! ! Conversion factor (ocean units are W/m2 and Kg/m2/s] zfacflx = 1.e0 ! no conversion [W/m2] ! W/m2 heat fluxes are send by the Atmosphere zfacwat = 1.e3 ! convert [m/s] to [kg/m2/s] ! m/s freshwater fluxes are send by the atmosphere ! ! =========================== ! ! ! Send Coupling fields ! ! ! =========================== ! ! !!gm bug ? here send instantaneous SST, not mean over the coupling period.... var_id = send_id(1) ; CALL cpl_prism_send( var_id, date, tn(:,:,1)+rt0, info ) ! ocean surface temperature [K] var_id = send_id(2) ; CALL cpl_prism_send( var_id, date, 1.0-frld , info ) ! fraction of ice-cover #if defined key_cpl_albedo DO jj = 1, jpj DO ji = 1, jpi IF( ( tn_ice(ji,jj) < 50 .OR. tn_ice(ji,jj) > 400 ) .AND. frld(ji,jj) < 1. ) THEN WRITE(numout,*) ' tn_ice, ERROR ', ji, jj, ' = ', tn_ice(ji,jj), & & ' qns_ice_recv=', qns_ice_recv(ji,jj), ' dqns_ice_recv=', dqns_ice_recv(ji,jj) ENDIF END DO END DO var_id = send_id(3) ; CALL cpl_prism_send( var_id, date, tn_ice , info ) ! ice surface temperature [K] var_id = send_id(4) ; CALL cpl_prism_send( var_id, date, alb_ice , info ) ! ice albedo [%] #else var_id = send_id(3) ; CALL cpl_prism_send( var_id, date, hicif , info ) ! ice thickness [m] var_id = send_id(4) ; CALL cpl_prism_send( var_id, date, hsnif , info ) ! snow thickness [m] #endif #if defined key_cpl_ocevel !!gm bug??? I have to check the grid point position... !! a priori there is a error here as un, vn are not at the same grid point.... !! there should be a averaged to set u and v at T-point.... with caution for sea-ice velocity at I-point.... un_weighted = un(:,:,1) * frld + u_ice * ( 1. - frld ) vn_weighted = vn(:,:,1) * frld + v_ice * ( 1. - frld ) ikchoix = - 1 ! converte from (i,j) to geographic referential CALL repere( un_weighted, vn_weighted, un_send, vn_send, ikchoix ) !!gm bug : at lbc_lnk is to be added on un_send and vn_send var_id = send_id(5) ; CALL cpl_prism_send( var_id, date, un_send , info ) ! surface current [m/s] var_id = send_id(6) ; CALL cpl_prism_send( var_id, date, vn_send , info ) ! surface current [m/s] #endif ! ! =========================== ! ! ! Recieve Momentum fluxes ! ! ! =========================== ! ! ! ... Receive wind stress fields in geographic component over water and ice var_id = recv_id(1) ; CALL cpl_prism_recv( var_id, date, ztaueuw, info ) ! ??? var_id = recv_id(2) ; CALL cpl_prism_recv( var_id, date, ztaunuw, info ) var_id = recv_id(3) ; CALL cpl_prism_recv( var_id, date, ztaueui, info ) var_id = recv_id(4) ; CALL cpl_prism_recv( var_id, date, ztaunui, info ) var_id = recv_id(5) ; CALL cpl_prism_recv( var_id, date, ztauevw, info ) var_id = recv_id(6) ; CALL cpl_prism_recv( var_id, date, ztaunvw, info ) var_id = recv_id(7) ; CALL cpl_prism_recv( var_id, date, ztauevi, info ) var_id = recv_id(8) ; CALL cpl_prism_recv( var_id, date, ztaunvi, info ) ! !!gm bug : keep separate ice and ocean stress ! ! ... combine water / ice stresses ztaueu(:,:) = ztaueuw(:,:) * frld(:,:) + ztaueui(:,:) * ( 1.0 - frld(:,:) ) ztaunu(:,:) = ztaunuw(:,:) * frld(:,:) + ztaunui(:,:) * ( 1.0 - frld(:,:) ) ztauev(:,:) = ztauevw(:,:) * frld(:,:) + ztauevi(:,:) * ( 1.0 - frld(:,:) ) ztaunv(:,:) = ztaunvw(:,:) * frld(:,:) + ztaunvi(:,:) * ( 1.0 - frld(:,:) ) ! ! ... rotate vector components from geographic to (i,j) referential CALL repcmo ( ztaueu, ztaunu, ztauev, ztaunv, utau, vtau, kt ) ! !!gm bug?? not sure but put that for security CALL lbc_lnk( utau , 'U', -1. ) CALL lbc_lnk( vtau , 'V', -1. ) !!gm end bug?? ! ! ! =========================== ! ! ! Recieve heat fluxes ! ! ! =========================== ! ! var_id = recv_id(13) ; CALL cpl_prism_recv( var_id, date, qsr_oce_recv , info ) ! ocean surface net downward shortwave flux var_id = recv_id(14) ; CALL cpl_prism_recv( var_id, date, qns_oce_recv , info ) ! ocean surface downward non-solar heat flux var_id = recv_id(15) ; CALL cpl_prism_recv( var_id, date, qsr_ice_recv , info ) ! ice solar heat flux var_id = recv_id(16) ; CALL cpl_prism_recv( var_id, date, qns_ice_recv , info ) ! ice non-solar heat flux var_id = recv_id(17) ; CALL cpl_prism_recv( var_id, date, dqns_ice_recv, info ) ! ice non-solar heat flux sensitivity qsr_oce_recv (:,:) = qsr_oce_recv (:,:) * tmask(:,:,1) * zfacflx qns_oce_recv (:,:) = qns_oce_recv (:,:) * tmask(:,:,1) * zfacflx qsr_ice_recv (:,:) = qsr_ice_recv (:,:) * tmask(:,:,1) * zfacflx qns_ice_recv (:,:) = qns_ice_recv (:,:) * tmask(:,:,1) * zfacflx dqns_ice_recv(:,:) = dqns_ice_recv(:,:) * tmask(:,:,1) * zfacflx IF( kt == nit000 ) THEN ! set once for all qsr penetration in sea-ice ! ! Since cloud cover catm not transmitted from atmosphere, it is set to 0. ! ! i.e. constant penetration fractions of 0.18 and 0.82 ! fraction of net shortwave radiation which is not absorbed in the thin surface layer and penetrates ! inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 ) fr1_i0_recv(:,:) = 0.18 fr2_i0_recv(:,:) = 0.82 ENDIF ! ! ! =========================== ! ! ! Recieve freshwater fluxes ! ! ! =========================== ! ! var_id = recv_id( 9) ; CALL cpl_prism_recv( var_id, date, zpew , info ) ! P-E over water var_id = recv_id(10) ; CALL cpl_prism_recv( var_id, date, zpei , info ) ! P-E over ice var_id = recv_id(11) ; CALL cpl_prism_recv( var_id, date, zpsol , info ) ! Snow fall over water and ice var_id = recv_id(12) ; CALL cpl_prism_recv( var_id, date, zevice, info ) ! Evaporation over ice (sublimination) ! ! ... calculate water flux (P-E over open ocean and ice) and solid precipitation (positive upward) tprecip_recv(:,:) = ( zpew (:,:) + zpei (:,:) ) * tmask(:,:,1) * zfacwat sprecip_recv(:,:) = ( zpsol(:,:) + zevice(:,:) ) * tmask(:,:,1) * zfacwat ! ... Control print & check IF(ln_ctl) THEN WRITE(numout,*) ' flx:tprecip_recv - Minimum value is ', MINVAL( tprecip_recv ) WRITE(numout,*) ' flx:tprecip_recv - Maximum value is ', MAXVAL( tprecip_recv ) WRITE(numout,*) ' flx:tprecip_recv - Sum value is ', SUM ( tprecip_recv ) ENDIF !!gm bug in mpp SUM require a mmp_sum call !!gm further more this test is quite expensive ... only needed at the first time-step??? IF( SUM( zpew*e1t*e2t ) /= SUM( zpew*e1t*e2t*tmask(:,:,1) ) ) THEN WRITE(numout,*) ' flx: Forcing values outside Orca mask' WRITE(numout,*) ' flx: Losses in water conservation' WRITE(numout,*) ' flx: Masked ', SUM(zpew*e1t*e2t*tmask(:,:,1)) WRITE(numout,*) ' flx: Unmasked ', SUM(zpew*e1t*e2t) WRITE(numout,*) ' flx: Simulation STOP' CALL FLUSH(numout) STOP END IF ! #if defined key_cpl_discharge ! Runoffs var_id = recv_id(18) ; CALL cpl_prism_recv ( var_id, date, calving_recv, info ) ! ice discharge into ocean var_id = recv_id(19) ; CALL cpl_prism_recv ( var_id, date, zrunriv , info ) ! river discharge into ocean var_id = recv_id(20) ; CALL cpl_prism_recv ( var_id, date, zruncot , info ) ! continental discharge into ocean DO jj = 1, jpj DO ji = 1, jpi zfact = zfacwat * tmask(ji,jj,1) calving_recv(ji,jj) = calving_recv(ji,jj) * zfact rrunoff_recv(ji,jj) = ( zrunriv(ji,jj) + zruncot(ji,jj) ) * zfact END DO END DO #else calving_recv(:,:) = 0. rrunoff_recv(:,:) = 0. #endif !!gm bug : this is not valid in mpp !!gm and I presum this is not required at all as a lbc_lnk is applied to all the fields at the end of the routine ! Oasis mask shift and update lateral boundary conditions (E. Maisonnave) ! not tested when mpp is used, W. Park !WSPTEST qsr_oce_recv (jpi-1,:) = qsr_oce_recv (1,:) qsr_ice_recv (jpi-1,:) = qsr_ice_recv (1,:) qns_oce_recv (jpi-1,:) = qns_oce_recv (1,:) qns_ice_recv (jpi-1,:) = qns_ice_recv (1,:) dqns_ice_recv(jpi-1,:) = dqns_ice_recv(1,:) tprecip_recv (jpi-1,:) = tprecip_recv (1,:) sprecip_recv (jpi-1,:) = sprecip_recv (1,:) fr1_i0_recv (jpi-1,:) = fr1_i0_recv (1,:) fr2_i0_recv (jpi-1,:) = fr2_i0_recv (1,:) rrunoff_recv (jpi-1,:) = rrunoff_recv (1,:) calving_recv (jpi-1,:) = calving_recv (1,:) !!gm end bug qsr (:,:) = qsr_oce_recv (:,:) ! ocean surface boundary condition qns (:,:) = qns_oce_recv (:,:) emp (:,:) = zpew (:,:) emps (:,:) = zpew (:,:) qsr_ice (:,:) = qsr_ice_recv (:,:) ! ice forcing fields qns_ice (:,:) = qns_ice_recv (:,:) dqns_ice(:,:) = dqns_ice_recv(:,:) tprecip (:,:) = tprecip_recv (:,:) sprecip (:,:) = sprecip_recv (:,:) fr1_i0 (:,:) = fr1_i0_recv (:,:) fr2_i0 (:,:) = fr2_i0_recv (:,:) !WSP rrunoff = rrunoff_recv !WSP calving = calving_recv rrunoff (:,:) = 0.e0 !WSP runoff and calving included in tprecip calving (:,:) = 0.e0 !WSP IF(ln_ctl) THEN WRITE(numout,*) 'flx:qsr_oce - Minimum value is ', MINVAL( qsr_oce ) WRITE(numout,*) 'flx:qsr_oce - Maximum value is ', MAXVAL( qsr_oce ) WRITE(numout,*) 'flx:qsr_oce - Sum value is ', SUM ( qsr_oce ) ! WRITE(numout,*) 'flx:tprecip - Minimum value is ', MINVAL( tprecip ) WRITE(numout,*) 'flx:tprecip - Maximum value is ', MAXVAL( tprecip ) WRITE(numout,*) 'flx:tprecip - Sum value is ', SUM ( tprecip ) ENDIF CALL lbc_lnk( qsr_oce , 'T', 1. ) CALL lbc_lnk( qsr_ice , 'T', 1. ) CALL lbc_lnk( qns_oce , 'T', 1. ) CALL lbc_lnk( qns_ice , 'T', 1. ) CALL lbc_lnk( tprecip , 'T', 1. ) CALL lbc_lnk( sprecip , 'T', 1. ) CALL lbc_lnk( rrunoff , 'T', 1. ) CALL lbc_lnk( dqns_ice, 'T', 1. ) CALL lbc_lnk( calving , 'T', 1. ) CALL lbc_lnk( fr1_i0 , 'T', 1. ) CALL lbc_lnk( fr2_i0 , 'T', 1. ) IF(ln_ctl) THEN WRITE(numout,*) 'flx(af lbc_lnk):qsr_oce - Minimum value is ', MINVAL( qsr_oce ) WRITE(numout,*) 'flx(af lbc_lnk):qsr_oce - Maximum value is ', MAXVAL( qsr_oce ) WRITE(numout,*) 'flx(af lbc_lnk):qsr_oce - Sum value is ', SUM ( qsr_oce ) ! WRITE(numout,*) 'flx(af lbc_lnk):tprecip - Minimum value is ', MINVAL( tprecip ) WRITE(numout,*) 'flx(af lbc_lnk):tprecip - Maximum value is ', MAXVAL( tprecip ) WRITE(numout,*) 'flx(af lbc_lnk):tprecip - Sum value is ', SUM ( tprecip ) ENDIF ! END SUBROUTINE sbc_cpl #else !!---------------------------------------------------------------------- !! Dummy routine NO sea surface restoring !!---------------------------------------------------------------------- LOGICAL, PUBLIC :: lk_sbc_cpl = .FALSE. !: coupled formulation flag CONTAINS SUBROUTINE sbc_cpl( kt ) ! Dummy routine WRITE(*,*) 'sbc_cpl: you should not have seen that print! error?', kt END SUBROUTINE sbc_cpl #endif !!====================================================================== END MODULE sbccpl