!Completed MODULE ocean_slab_mod ! !================================================================== ! ! Purpose ! ------- ! The dynamical slab ocean model of the Generic-PCM. It has the following features: ! (a) Computes sea ice creation and evolution. ! (b) Snow has thermodynamic properties. ! (c) Computes oceanic horizontal transport (diffusion & surface-wind driven Ekman transport). ! (d) Can be used in parallel mode. ! ! Authors ! ------- ! S. Bhatnagar and E. Millour (2023) ! Adapted from the ocean modules of LMDZ Earth (F. Codron) and the Generic-PCM (B. Charnay, 2013). ! ! Notes ! ----- ! Compared to the old model, the new model has the following changes (non-exhaustive): ! (a) More realistic description of sea ice creation and evolution - simultaneous ! surface, side and bottom melting / freezing depending on fluxes. ! (b) Snow has an effective heat capacity. ! (c) Snow has "weight"; it can sink an ice block if there is too much of it. ! (d) Snow can be blown off by wind. ! (e) The two-layer ocean allows for convective adjustment. ! (f) Diffusion can follow the Gent-McWilliams scheme + Eddy diffusivity. ! (g) Can be used in parallel mode. ! !================================================================== USE dimphy, ONLY: klon USE mod_grid_phy_lmdz, ONLY: klon_glo USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root IMPLICIT NONE PRIVATE PUBLIC :: ocean_slab_init, ocean_slab_ice, ocean_slab_noice, & ocean_slab_frac, ocean_slab_get_vars, ocean_slab_final !*********************************************************************************** ! Global saved variables !*********************************************************************************** ! number of slab vertical layers INTEGER, PUBLIC, SAVE :: nslay=2 !$OMP THREADPRIVATE(nslay) ! number of oceanic grid points INTEGER, PUBLIC, SAVE :: knon !$OMP THREADPRIVATE(knon) ! timestep for coupling (update slab temperature) in timesteps INTEGER, PRIVATE, SAVE :: cpl_pas !$OMP THREADPRIVATE(cpl_pas) ! cyang = 1/heat capacity of top layer (rho.c.H) REAL, PRIVATE, SAVE :: cyang !$OMP THREADPRIVATE(cyang) ! depth of slab layers (1st or 2nd layer) REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: slabh !$OMP THREADPRIVATE(slabh) ! slab temperature REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: tslab !$OMP THREADPRIVATE(tslab) ! heat flux convergence due to Ekman REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_ekman !$OMP THREADPRIVATE(dt_ekman) ! heat flux convergence due to horiz diffusion REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_hdiff !$OMP THREADPRIVATE(dt_hdiff) ! heat flux convergence due to GM eddy advection REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_gm !$OMP THREADPRIVATE(dt_gm) ! fraction of ocean covered by sea ice (sic / (oce+sic)) REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: fsic !$OMP THREADPRIVATE(fsic) ! temperature of the sea ice REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: tice !$OMP THREADPRIVATE(tice) ! sea ice thickness, in kg/m2 REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: seaice !$OMP THREADPRIVATE(seaice) ! net surface heat flux, weighted by open ocean fraction ! slab_bils accumulated over cpl_pas timesteps REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bils_cum !$OMP THREADPRIVATE(bils_cum) ! net heat flux into the ocean below the ice : conduction + solar radiation REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: slab_bilg !$OMP THREADPRIVATE(slab_bilg) ! slab_bilg cululated over cpl_pas timesteps REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bilg_cum !$OMP THREADPRIVATE(bilg_cum) ! wind stress saved over cpl_pas timesteps REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: taux_cum !$OMP THREADPRIVATE(taux_cum) REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: tauy_cum !$OMP THREADPRIVATE(tauy_cum) !*********************************************************************************** ! Parameters (could be read in def file: move to slab_init) !*********************************************************************************** ! snow and ice physical characteristics: REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp [in K] REAL, PARAMETER :: t_melt=273.15 ! melting ice temp [in K] REAL, PARAMETER :: sno_den=300. !mean snow density [in kg/m3] REAL, PARAMETER :: ice_den=917. ! ice density [in kg/m3] REAL, PARAMETER :: sea_den=1026. ! sea water density [in kg/m3] REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity of ice [in W/(m.K) or (W.kg)/(K.m4)] REAL, PRIVATE, SAVE :: sno_cond ! conductivity of snow [in W/(m.K) or (W.kg)/(K.m4)] !$OMP THREADPRIVATE(sno_cond) REAL, PARAMETER :: ice_cap=2067. ! specific heat capacity, snow and ice [in J/(kg.K)] REAL, PARAMETER :: sea_cap=3994. ! specific heat capacity, seawater [in J/(kg.K)] REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice [in J/kg] REAL, PARAMETER :: ice_sub=2834000. ! latent heat of sublimation for snow and ice [in J/kg] ! control of snow and ice cover & freeze / melt (heights in m converted to kg/m2) REAL, PARAMETER, PUBLIC :: snow_min=0.05*sno_den ! critical snow height [in kg/m2] REAL, PARAMETER :: snow_wfact=0.4 ! max fraction of falling snow blown into ocean [in kg/m2] REAL, PARAMETER :: ice_frac_min=0.001 REAL, PRIVATE, SAVE :: ice_frac_max ! Max ice fraction (leads) !$OMP THREADPRIVATE(ice_frac_max) REAL, PARAMETER :: h_ice_min=0.01*ice_den ! min ice thickness [in kg/m2] REAL, PARAMETER :: h_ice_thin=0.15*ice_den ! thin ice thickness [in kg/m2] ! below ice_thin, priority is to melt lateral / grow height ! ice_thin is also height of new ice REAL, PRIVATE, SAVE :: h_ice_thick ! thin ice thickness !$OMP THREADPRIVATE(h_ice_thick) ! above ice_thick, priority is melt height / grow lateral REAL, PARAMETER :: h_ice_new=1.*ice_den ! max height of new open ocean ice [in kg/m2] REAL, PARAMETER :: h_ice_max=10.*ice_den ! max ice height [in kg/m2] REAL, PARAMETER :: epsfra=1.0E-05 ! minimial grid fraction size below which there is no ice REAL, PARAMETER, PUBLIC :: capcalocean=50.*4.228e+06 ! surface heat capacity [J.K-1.m-2] (assuming 50 m slab ocean) REAL, PARAMETER, PUBLIC :: capcalseaice=5.1444e+06*0.15 REAL, PARAMETER, PUBLIC :: capcalsno=2.3867e+06*0.15 REAL, PARAMETER, PUBLIC :: h_alb_ice=0.3*ice_den ! height (in kg/m2) used in the calculation of sea ice albedo vs thickness ! (changed from 50cm to 30cm based on comparisons with Brandt et al. 2005) ; more info in the slab ocean wiki page REAL, PARAMETER, PUBLIC :: h_sno_alb=0.02*sno_den ! height (in kg/m2) for control of snow fraction REAL, PARAMETER, PUBLIC :: alb_ice_min=0.08 ! minimum sea ice albedo used for calculation of albedo as a function of sea ice thickness (https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Slab_ocean_model) ! Horizontal transport parameters ! flag for horizontal diffusion LOGICAL, PUBLIC, SAVE :: slab_hdiff !$OMP THREADPRIVATE(slab_hdiff) ! flag for GM eddy diffusivity LOGICAL, PUBLIC, SAVE :: slab_gm !$OMP THREADPRIVATE(slab_gm) REAL, PRIVATE, SAVE :: coef_hdiff ! coefficient for horizontal diffusion !$OMP THREADPRIVATE(coef_hdiff) ! flags for Ekman, conv adjustment LOGICAL, PUBLIC, SAVE :: slab_ekman !$OMP THREADPRIVATE(slab_ekman) INTEGER, PUBLIC, SAVE :: slab_cadj !$OMP THREADPRIVATE(slab_cadj) !*********************************************************************************** CONTAINS ! !*********************************************************************************** ! SUBROUTINE ocean_slab_init(dtime, pctsrf_rst, tslab_rst, tice_rst, seaice_rst, zmasq) ! This routine ! (1) allocates variables initialised from restart fields ! (2) allocates some other variables internal to the ocean module USE ioipsl_getin_p_mod, ONLY : getin_p USE mod_phys_lmdz_transfert_para, ONLY : gather USE slab_heat_transp_mod, ONLY : ini_slab_transp ! Input variables !*********************************************************************************** REAL, INTENT(IN) :: dtime ! Variables read from restart file REAL, DIMENSION(klon), INTENT(IN) :: pctsrf_rst REAL, DIMENSION(klon,nslay), INTENT(IN) :: tslab_rst REAL, DIMENSION(klon), INTENT(IN) :: tice_rst REAL, DIMENSION(klon), INTENT(IN) :: seaice_rst REAL, DIMENSION(klon), INTENT(IN) :: zmasq ! Local variables !************************************************************************************ INTEGER :: error REAL, DIMENSION(klon_glo) :: zmasq_glo CHARACTER (len = 80) :: abort_message CHARACTER (len = 20) :: modname = 'ocean_slab_init' !*********************************************************************************** ! Define some parameters !*********************************************************************************** ! ! cpl_pas coupling period (update of tslab and ice fraction) ! for a calculation at each physical timestep, cpl_pas=1 cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour CALL getin_p('cpl_pas',cpl_pas) print *,'cpl_pas',cpl_pas ! Number of slab layers ! nslay=2 ! CALL getin_p('slab_layers',nslay) print *,'number of slab layers : ',nslay ! Layer thickness ALLOCATE(slabh(nslay), stat = error) IF (error /= 0) THEN abort_message='Pb allocation slabh' CALL abort_physic(modname,abort_message,1) ENDIF slabh(1)=50. ! Height of first ocean layer (wind-mixed layer) CALL getin_p('slab_depth',slabh(1)) IF (nslay.GT.1) THEN slabh(2)=150. ! Height of second ocean layer (deep ocean layer) END IF ! cyang = 1/heat capacity of top layer (rho.c.H) cyang=1/(slabh(1)*sea_den*sea_cap) ! ********** Sea Ice parameters *********** ice_frac_max = 0.999999 ! frac = 1 may lead to some problems. CALL getin_p('ice_frac_max',ice_frac_max) h_ice_thick = 1.5 CALL getin_p('h_ice_thick',h_ice_thick) h_ice_thick = h_ice_thick * ice_den sno_cond = 0.31 CALL getin_p('sno_cond',sno_cond) sno_cond = sno_cond * sno_den ! ********** Heat Transport parameters **** ! Ekman transport ! slab_ekman=0 slab_ekman=.FALSE. CALL getin_p('slab_ekman',slab_ekman) ! GM eddy advection (2-layers only) slab_gm=.FALSE. CALL getin_p('slab_gm',slab_gm) ! IF (slab_ekman.LT.2) THEN IF (.NOT.slab_ekman) THEN slab_gm=.FALSE. ENDIF ! Horizontal diffusion slab_hdiff=.FALSE. CALL getin_p('slab_hdiff',slab_hdiff) IF (slab_gm) THEN coef_hdiff=8000. ! non-dimensional; coef_hdiff should be 25000 if GM is off ELSE coef_hdiff=25000. ENDIF CALL getin_p('coef_hdiff',coef_hdiff) ! Convective adjustment ! IF (nslay.EQ.1) THEN ! slab_cadj=0 ! ELSE slab_cadj=1 ! END IF CALL getin_p('slab_cadj',slab_cadj) !************************************************************************************ ! Allocate surface fraction read from restart file !************************************************************************************ ALLOCATE(fsic(klon), stat = error) IF (error /= 0) THEN abort_message='Pb allocation tmp_pctsrf_slab' CALL abort_physic(modname,abort_message,1) ENDIF fsic(:)=0. !zmasq = continent fraction WHERE (1.-zmasq(:)>EPSFRA) ! fsic(:) = MIN(pctsrf_rst(:,is_sic)/(1.-zmasq(:)),ice_frac_max) fsic(:) = MIN(pctsrf_rst(:)/(1.-zmasq(:)),ice_frac_max) END WHERE !************************************************************************************ ! Allocate saved fields !************************************************************************************ ALLOCATE(tslab(klon,nslay), stat=error) IF (error /= 0) CALL abort_physic & (modname,'pb allocation tslab', 1) tslab(:,:) = tslab_rst(:,:) ALLOCATE(bils_cum(klon), stat = error) IF (error /= 0) THEN abort_message='Pb allocation slab_bils_cum' CALL abort_physic(modname,abort_message,1) ENDIF bils_cum(:) = 0.0 ! IF (version_ocean=='sicINT') THEN ! interactive sea ice ALLOCATE(slab_bilg(klon), stat = error) IF (error /= 0) THEN abort_message='Pb allocation slab_bilg' CALL abort_physic(modname,abort_message,1) ENDIF slab_bilg(:) = 0.0 ALLOCATE(bilg_cum(klon), stat = error) IF (error /= 0) THEN abort_message='Pb allocation slab_bilg_cum' CALL abort_physic(modname,abort_message,1) ENDIF bilg_cum(:) = 0.0 ALLOCATE(tice(klon), stat = error) IF (error /= 0) THEN abort_message='Pb allocation slab_tice' CALL abort_physic(modname,abort_message,1) ENDIF tice(:) = tice_rst(:) ALLOCATE(seaice(klon), stat = error) IF (error /= 0) THEN abort_message='Pb allocation slab_seaice' CALL abort_physic(modname,abort_message,1) ENDIF seaice(:) = seaice_rst(:) ! END IF IF (slab_hdiff) THEN !horizontal diffusion ALLOCATE(dt_hdiff(klon,nslay), stat = error) IF (error /= 0) THEN abort_message='Pb allocation dt_hdiff' CALL abort_physic(modname,abort_message,1) ENDIF dt_hdiff(:,:) = 0.0 ENDIF IF (slab_gm) THEN !GM advection ALLOCATE(dt_gm(klon,nslay), stat = error) IF (error /= 0) THEN abort_message='Pb allocation dt_gm' CALL abort_physic(modname,abort_message,1) ENDIF dt_gm(:,:) = 0.0 ENDIF ! IF (slab_ekman.GT.0) THEN ! ekman transport IF (slab_ekman) THEN ! ekman transport ALLOCATE(dt_ekman(klon,nslay), stat = error) IF (error /= 0) THEN abort_message='Pb allocation dt_ekman' CALL abort_physic(modname,abort_message,1) ENDIF dt_ekman(:,:) = 0.0 ALLOCATE(taux_cum(klon), stat = error) IF (error /= 0) THEN abort_message='Pb allocation taux_cum' CALL abort_physic(modname,abort_message,1) ENDIF taux_cum(:) = 0.0 ALLOCATE(tauy_cum(klon), stat = error) IF (error /= 0) THEN abort_message='Pb allocation tauy_cum' CALL abort_physic(modname,abort_message,1) ENDIF tauy_cum(:) = 0.0 ENDIF ! Initialize transport IF (slab_hdiff.OR.(slab_ekman)) THEN CALL gather(zmasq,zmasq_glo) ! Master thread/process only !$OMP MASTER IF (is_mpi_root) THEN CALL ini_slab_transp(zmasq_glo) END IF !$OMP END MASTER END IF END SUBROUTINE ocean_slab_init ! !*********************************************************************************** ! SUBROUTINE ocean_slab_frac(pctsrf_chg, zmasq) ! This routine sends back the sea ice and ocean fraction to the main physics routine. ! Called only with interactive sea ice. ! Arguments !************************************************************************************ REAL, DIMENSION(klon), INTENT(IN) :: zmasq ! proxy of rnat REAL, DIMENSION(klon), INTENT(OUT) :: pctsrf_chg ! sea-ice fraction pctsrf_chg(:)=fsic(:)*(1.-zmasq(:)) END SUBROUTINE ocean_slab_frac ! !************************************************************************************ ! SUBROUTINE ocean_slab_noice(itime, dtime, knon, knindex, & precip_snow, tsurf_in, & radsol, snow, fluxsens, & tsurf_new, flux_u1, flux_v1, zmasq) USE slab_heat_transp_mod, ONLY: divgrad_phy,slab_ekman2,slab_gmdiff USE mod_phys_lmdz_para ! This routine ! (1) computes surface turbulent fluxes over points with some open ocean ! (2) reads additional Q-flux (everywhere) ! (3) computes horizontal transport (diffusion & Ekman) ! (4) updates slab temperature every cpl_pas ; creates new ice if needed. ! Note : ! klon total number of points ! knon number of points with open ocean (varies with time) ! knindex gives position of the knon points within klon. ! In general, local saved variables have klon values ! variables exchanged with PBL module have knon. ! Input arguments !*********************************************************************************** INTEGER, INTENT(IN) :: itime ! current timestep INTEGER, INTEGER, INTENT(IN) :: knon ! number of points INTEGER, DIMENSION(klon), INTENT(IN) :: knindex REAL, INTENT(IN) :: dtime ! timestep (s) REAL, DIMENSION(klon), INTENT(IN) :: precip_snow !, precip_rain REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in ! surface temperature REAL, DIMENSION(klon), INTENT(IN) :: radsol ! net surface (radiative) flux REAL, DIMENSION(klon), INTENT(IN) :: flux_u1, flux_v1 ! Comes from turbdiff REAL, DIMENSION(klon), INTENT(IN) :: fluxsens !, sensible heat flux REAL, DIMENSION(klon), INTENT(IN) :: zmasq ! proxy of rnat ! In/Output arguments !************************************************************************************ REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! in kg/m2 ! Output arguments !************************************************************************************ REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! new surface tempearture ! Local variables !************************************************************************************ INTEGER :: i,ki,k REAL :: t_cadj ! for new ice creation REAL :: e_freeze REAL, DIMENSION(klon) :: slab_bils ! weighted surface heat flux ! horizontal diffusion and Ekman local vars ! dimension = global domain (klon_glo) instead of // subdomains REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo,dt_ekman_glo,dt_gm_glo ! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_tmp, dt_ekman_tmp REAL, DIMENSION(klon_glo,nslay) :: tslab_glo REAL, DIMENSION(klon_glo) :: taux_glo,tauy_glo !**************************************************************************************** ! 1) Surface fluxes calculation ! Points with some open ocean only !**************************************************************************************** ! save total cumulated heat fluxes locally ! radiative + turbulent + melt of falling snow slab_bils(:)=0. DO i=1,knon ki=knindex(i) slab_bils(ki)=(1.-fsic(ki))*(fluxsens(ki)+radsol(ki) & -precip_snow(ki)*ice_lat*(1.+snow_wfact*fsic(ki))) bils_cum(ki)=bils_cum(ki)+slab_bils(ki) END DO ! Compute surface wind stress ! CALL calcul_flux_wind(knon, dtime, & ! u0, v0, u1, v1, gustiness, cdragm, & ! flux_u1, flux_v1) ! save cumulated wind stress IF (slab_ekman) THEN DO i=1,knon ki=knindex(i) taux_cum(ki)=taux_cum(ki)+flux_u1(ki)*(1.-fsic(ki))/cpl_pas tauy_cum(ki)=tauy_cum(ki)+flux_v1(ki)*(1.-fsic(ki))/cpl_pas END DO ENDIF !**************************************************************************************** ! 2) Q-Flux : get global variables lmt_bils, diff_sst and diff_siv from file ! limit_slab.nc ! !**************************************************************************************** ! CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) ! lmt_bils and diff_sst,siv saved by limit_slab ! qflux = total QFlux correction (in W/m2) ! IF (qflux_bnd.EQ.2) THEN ! dt_qflux(:,1) = lmt_bils(:,1)+diff_sst(:)/cyang/86400. ! dt_qflux_sic(:) = -diff_siv(:)*ice_den*ice_lat/86400. ! ELSE ! dt_qflux(:,1) = lmt_bils(:,1)+diff_sst(:)/cyang/86400. & ! - diff_siv(:)*ice_den*ice_lat/86400. ! END IF ! IF (nslay.GT.1) THEN ! dt_qflux(:,2:nslay)=lmt_bils(:,2:nslay) ! END IF !**************************************************************************************** ! 3) Recalculate new temperature (add Surf fluxes, Q-Flux, Ocean transport) ! Bring to freezing temp and make sea ice if necessary ! !***********************************************o***************************************** tsurf_new=tsurf_in IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction ! *********************************** ! Horizontal transport ! *********************************** IF (slab_ekman) THEN ! copy wind stress to global var CALL gather(taux_cum,taux_glo) CALL gather(tauy_cum,tauy_glo) END IF IF (slab_hdiff.OR.(slab_ekman)) THEN CALL gather(tslab,tslab_glo) ! Compute horiz transport on one process only IF (is_mpi_root .AND. is_omp_root) THEN ! Only master processus IF (slab_hdiff) THEN dt_hdiff_glo(:,:)=0. END IF IF (slab_ekman) THEN dt_ekman_glo(:,:)=0. END IF IF (slab_gm) THEN dt_gm_glo(:,:)=0. END IF DO i=1,cpl_pas ! time splitting for numerical stability ! IF (slab_ekman.GT.0) THEN ! SELECT CASE (slab_ekman) ! CASE (1) ! 1.5 layer scheme ! CALL slab_ekman1(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp) ! CASE (2) ! 2-layers ! CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp,dt_hdiff_tmp,slab_gm) ! CASE DEFAULT ! dt_ekman_tmp(:,:)=0. ! END SELECT IF (slab_ekman) THEN CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp,dt_hdiff_tmp,slab_gm) dt_ekman_glo(:,:)=dt_ekman_glo(:,:)+dt_ekman_tmp(:,:) ! convert dt_ekman from K.s-1.(kg.m-2) to K.s-1 DO k=1,nslay dt_ekman_tmp(:,k)=dt_ekman_tmp(:,k)/(slabh(k)*sea_den) ENDDO tslab_glo=tslab_glo+dt_ekman_tmp*dtime IF (slab_gm) THEN ! Gent-McWilliams eddy advection dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:) ! convert dt from K.s-1.(kg.m-2) to K.s-1 DO k=1,nslay dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/(slabh(k)*sea_den) END DO tslab_glo=tslab_glo+dt_hdiff_tmp*dtime END IF ENDIF IF (slab_hdiff) THEN ! horizontal diffusion ! laplacian of slab T CALL divgrad_phy(nslay,tslab_glo,dt_hdiff_tmp) ! multiply by diff coef and normalize to 50m slab equivalent dt_hdiff_tmp=dt_hdiff_tmp*coef_hdiff*50./SUM(slabh) dt_hdiff_glo(:,:)=dt_hdiff_glo(:,:)+ dt_hdiff_tmp(:,:) tslab_glo=tslab_glo+dt_hdiff_tmp*dtime END IF END DO ! time splitting IF (slab_hdiff) THEN !dt_hdiff_glo saved in W/m2 DO k=1,nslay dt_hdiff_glo(:,k)=dt_hdiff_glo(:,k)*slabh(k)*sea_den*sea_cap/cpl_pas END DO END IF IF (slab_gm) THEN !dt_hdiff_glo saved in W/m2 dt_gm_glo(:,:)=dt_gm_glo(:,:)*sea_cap/cpl_pas END IF IF (slab_ekman) THEN ! dt_ekman_glo saved in W/m2 dt_ekman_glo(:,:)=dt_ekman_glo(:,:)*sea_cap/cpl_pas END IF END IF ! master process !$OMP BARRIER ! Send new fields back to all processes CALL Scatter(tslab_glo,tslab) IF (slab_hdiff) THEN CALL Scatter(dt_hdiff_glo,dt_hdiff) END IF IF (slab_gm) THEN CALL Scatter(dt_gm_glo,dt_gm) END IF IF (slab_ekman) THEN CALL Scatter(dt_ekman_glo,dt_ekman) ! clear wind stress taux_cum(:)=0. tauy_cum(:)=0. END IF ENDIF ! transport ! *********************************** ! Other heat fluxes ! *********************************** ! Add cumulated ocean surface fluxes tslab(:,1) = tslab(:,1) + bils_cum(:)*cyang*dtime ! Convective adjustment if 2 layers IF ((nslay.GT.1).AND.(slab_cadj.GT.0)) THEN DO i=1,klon IF (tslab(i,2).GT.tslab(i,1)) THEN ! mean (mass-weighted) temperature t_cadj=SUM(tslab(i,:)*slabh(:))/SUM(slabh(:)) tslab(i,1)=t_cadj tslab(i,2)=t_cadj END IF END DO END IF ! Add read QFlux ! IF (qflux_bnd.EQ.1) THEN ! ! QFlux from ocean circ. cannot cool tslab below freezing. ! dq_freeze = (t_freeze - tslab(:,1)) / (cyang*dtime*cpl_pas) ! dt_qflux(:,1) = MAX(dt_qflux(:,1), MIN(dq_freeze,0.)) ! ELSE IF (qflux_bnd.EQ.2) THEN ! dq_freeze = (t_freeze - tslab(:,1)) / (cyang*dtime*cpl_pas) & ! - dt_qflux_sic(:) ! dt_qflux(:,1) = MAX(dt_qflux(:,1), MIN(dq_freeze,0.)) & ! + dt_qflux_sic(:) ! END IF ! DO k=1,nslay ! tslab(:,k) = tslab(:,k) + dt_qflux(:,k)*cyang*dtime*cpl_pas & ! * slabh(1)/slabh(k) ! END DO ! *********************************** ! Update surface temperature and ice ! *********************************** ! SELECT CASE(version_ocean) ! CASE('sicNO') ! no sea ice even below freezing ! ! DO i=1,knon ! ki=knindex(i) ! tsurf_new(i)=tslab(ki,1) ! END DO ! CASE('sicOBS') ! "realistic" case, for prescribed sea ice ! ! tslab cannot be below freezing ! DO i=1,knon ! ki=knindex(i) ! IF (tslab(ki,1).LT.t_freeze) THEN ! tslab(ki,1)=t_freeze ! END IF ! tsurf_new(i)=tslab(ki,1) ! END DO ! CASE('sicINT') ! interactive sea ice DO i=1,knon ki=knindex(i) ! Check if new slab temperature is below freezing IF (tslab(ki,1).LT.t_freeze) THEN ! We need to form ice now over ice-free points ! Else points not seen by slab_ice IF (fsic(ki)*(1.-zmasq(ki)).LT.epsfra) THEN ! No ice present yet. ! quantity of new ice formed e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat & +fsic(ki)*seaice(ki) ! new ice forms at h_ice_thin tsurf_new(ki)=t_freeze tice(ki)=t_freeze fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin) IF (fsic(ki).GT.ice_frac_min) THEN seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max) tslab(ki,1)=t_freeze ELSE fsic(ki)=0. END IF END IF ! sea ice present ! if ice already present, new ice formed in slab_ice routine. ! ELSE ! temperature above freezing ! tsurf_new(i)=tslab(ki,1) END IF END DO ! END SELECT bils_cum(:)=0.0! clear cumulated fluxes END IF ! coupling time END SUBROUTINE ocean_slab_noice ! !***************************************************************************** ! SUBROUTINE ocean_slab_ice( & ! itime, dtime, jour, knon, knindex, & ! tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, & ! AcoefH, AcoefQ, BcoefH, BcoefQ, & ! AcoefU, AcoefV, BcoefU, BcoefV, & ! ps, u1, v1, gustiness, & ! radsol, snow, qsurf, qsol, agesno, & ! alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, & ! tsurf_new, dflux_s, dflux_l, swnet) SUBROUTINE ocean_slab_ice(itime, dtime, knon, knindex, & precip_snow, tsurf_in, & radsol, snow, fluxsens, & tsurf_new, evap, flux_u1, flux_v1, zmasq) ! USE calcul_fluxs_mod ! INCLUDE "YOMCST.h" ! INCLUDE "clesphys.h" ! This routine ! (1) computes surface turbulent fluxes over points with some sea ice ! (2) computes condutive fluxes in the snow and ice, and ice-ocean flux ! (3) computes snow/ice albedo ! (4) updates snow/ice temperature, surface melt if needed. ! (5) lateral growth if tslab < freezing ! (6) bottom & side melt / growth depending on bottom fluxes ! (7) updates slab temperature every cpl_pas ! Note : ! klon total number of points ! knon number of points with open ocean (varies with time) ! knindex gives position of the knon points within klon. ! In general, local saved variables have klon values ! variables exchanged with PBL module have knon. ! Input arguments !**************************************************************************************** INTEGER, INTENT(IN) :: itime, knon !, jour INTEGER, DIMENSION(klon), INTENT(IN) :: knindex REAL, INTENT(IN) :: dtime REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in ! REAL, DIMENSION(klon), INTENT(IN) :: p1lay ! REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm REAL, DIMENSION(klon), INTENT(IN) :: precip_snow !, precip_rain REAL, DIMENSION(klon), INTENT(IN) :: evap, fluxsens!,fluxlat REAL, DIMENSION(klon), INTENT(IN) :: flux_u1, flux_v1 REAL, DIMENSION(klon), INTENT(IN) :: zmasq ! proxy of rnat ! REAL, DIMENSION(klon), INTENT(IN) :: spechum, temp_air ! REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ ! REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV ! REAL, DIMENSION(klon), INTENT(IN) :: ps ! REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness ! REAL, DIMENSION(klon), INTENT(IN) :: swnet ! In/Output arguments !**************************************************************************************** REAL, DIMENSION(klon), INTENT(INOUT) :: snow!, qsol ! REAL, DIMENSION(klon), INTENT(INOUT) :: agesno REAL, DIMENSION(klon), INTENT(INOUT) :: radsol ! Output arguments !**************************************************************************************** ! REAL, DIMENSION(klon), INTENT(OUT) :: qsurf ! REAL, DIMENSION(klon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval ! REAL, DIMENSION(klon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval ! REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens!, fluxlat ! REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1 REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l ! Local variables !**************************************************************************************** INTEGER :: i,ki ! REAL, DIMENSION(klon) :: cal, beta, dif_grnd ! REAL, DIMENSION(klon) :: u0, v0 ! REAL, DIMENSION(klon) :: u1_lay, v1_lay REAL, DIMENSION(klon) :: f_bot ! bottom ocean - ice flux ! intermediate heat fluxes: ! (conduction snow / ice, shortwave penetration, bottom turbulent fluxes) REAL :: f_cond_s, f_cond_i, f_turb ! friction velocity, 1/C and 1/tau conduction for ice REAL :: ustar ! REAL :: uscap, ustau ! for snow/ice albedo: ! REAL :: alb_snow, alb_ice, alb_pond ! REAL :: frac_snow, frac_ice, frac_pond ! REAL :: z1_i, z2_i, z1_s, zlog ! height parameters ! for ice melt / freeze REAL :: e_melt, e_freeze, snow_evap, h_test, h_new ! dhsic, dfsic change in ice mass, fraction. ! frac_mf ratio of lateral / thickness growth / melt REAL :: dhsic, dfsic, frac_mf !******************************************************************************* ! 1) Update surface , ice and slab temperature !******************************************************************************* ! Wind stress ! flux_u1, flux_v1 from physics ! save cumulated wind stress ! Use ocean-ice stress = wind stress * (1.-fsic) ! IF (slab_ekman.GT.0) THEN IF (slab_ekman) THEN DO i=1,knon ki=knindex(i) taux_cum(ki)=taux_cum(ki)+flux_u1(ki)*fsic(ki)*(1.-fsic(ki))/cpl_pas tauy_cum(ki)=tauy_cum(ki)+flux_v1(ki)*fsic(ki)*(1.-fsic(ki))/cpl_pas END DO ENDIF ! Initialize ice-ocean flux slab_bilg(:)=0. ! Old, explicit scheme for snow & ice conductive fluxes ! radsol is total surface fluxes (input) : radiative + turbulent DO i=1,knon ki=knindex(i) ! For PCM : you can probably replace ki by i ! ocean-ice turbulent heat flux ! turbulent velocity = (tau/rho)^1/2 ustar = MAX(5e-4, SQRT((1.-fsic(ki)) & * SQRT(flux_u1(ki)**2 + flux_v1(ki)**2) / sea_den)) f_turb = 0.0057 * sea_den * sea_cap * ustar * (tslab(ki,1) - t_freeze) ! f_turb >0 and cannot bring tslab below zero f_turb = MAX(0., MIN(f_turb, & (tslab(ki,1)-t_freeze) / (cyang*dtime*cpl_pas))) ! ice conductive flux (pos up) IF (seaice(ki).GT.0) THEN f_cond_i = ice_cond*(t_freeze-tice(ki))/seaice(ki) ELSE f_cond_i = 0 END IF ! If snow layer present, tsurf = tsnow ! Problem here, as tsurf_in # tsnow ? IF (snow(ki).GT.snow_min) THEN ! snow conductive flux (pos up) f_cond_s=sno_cond*(tice(ki)-tsurf_in(ki))/snow(ki) ! update ice temperature tice(ki)=tice(ki) + 2./ice_cap/(snow(ki)+seaice(ki)) & *(f_cond_i-f_cond_s)*dtime ! update snow temperature tsurf_new(ki) = tsurf_in(ki) + 2./ice_cap/snow(ki) & *(fluxsens(ki)+radsol(ki)+f_cond_s)*dtime ELSE IF (seaice(ki).GT.0) THEN ! bare ice. tsurf = tice ! update ice temperature tice(ki) = tice(ki) + 2./ice_cap/seaice(ki) & *(fluxsens(ki)+radsol(ki)+f_cond_i)*dtime tsurf_new(ki) = tice(ki) END IF ! bottom flux (used to grow ice from below) f_bot(ki) = f_turb - f_cond_i slab_bilg(ki) = -f_turb END DO ! !! Surface turbulent fluxes (sens, lat etc) and update surface temp. ! dif_grnd(:)=0. ! beta(:) = 1. ! CALL calcul_fluxs(knon, is_sic, dtime, & ! tsurf_in, p1lay, cal, beta, cdragh, cdragh, ps, & ! precip_rain, precip_snow, snow, qsurf, & ! radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & ! f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, & ! tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l) ! DO i=1,knon ! IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i) ! END DO ! Surface, snow-ice and ice-ocean fluxes. ! Prepare call to calcul_fluxs (cal, beta, radsol, dif_grnd) ! Surface turbulent fluxes (sens, lat etc) and update surface temp. ! beta(:) = 1. ! CALL calcul_fluxs(knon, is_sic, dtime, & ! tsurf_in, p1lay, cal, beta, cdragh, cdragh, ps, & ! precip_rain, precip_snow, snow, qsurf, & ! radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & ! f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, & ! tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l) !! Update remaining temperature and fluxes ! DO i=1,knon ! ki=knindex(i) ! ! ocean-ice turbulent heat flux ! ! turbulent velocity = (tau/rho)^1/2 for low ice cover ! ! min = 5e-4 for 1cm/s current ! ustar = MAX(5e-4, SQRT((1.-fsic(ki)) & ! * SQRT(flux_u1(i)**2 + flux_v1(i)**2) / sea_den)) ! f_turb = 0.0057 * sea_den * sea_cap * ustar * (tslab(ki,1) - t_freeze) ! ! ice_ocean fluxes cannot bring tslab below freezing ! f_turb = MAX(0., MIN(f_turb, slab_bilg(ki) + (tslab(ki,1)-t_freeze) & ! / (fsic(ki)*cyang*dtime*cpl_pas))) ! IF (snow(i).GT.snow_min) THEN ! ! snow conductive flux after calcul_fluxs ! f_cond_s = sno_cond * (tice(ki)-tsurf_new(i)) / snow(i) ! ! 1 / heat capacity and conductive timescale ! uscap = 2. / ice_cap / (snow(i)+seaice(ki)) ! ustau = uscap * ice_cond / seaice(ki) ! ! update ice temp ! tice(ki) = (tice(ki) + dtime*(ustau*t_freeze - uscap*f_cond_s)) & ! / (1 + dtime*ustau) ! ELSE ! bare ice ! tice(ki)=tsurf_new(i) ! END IF ! ! ice conductive flux (pos up) ! f_cond_i = ice_cond * (t_freeze-tice(ki)) / seaice(ki) ! f_bot(i) = f_turb - f_cond_i ! slab_bilg(ki) = slab_bilg(ki)-f_turb ! END DO ! weight fluxes to ocean by sea ice fraction slab_bilg(:)=slab_bilg(:)*fsic(:) !**************************************************************************************** ! 2) Update snow and ice surface : thickness and fraction !**************************************************************************************** DO i=1,knon ki=knindex(i) ! snow precip (could be before fluxes above ?) IF (precip_snow(ki) > 0.) THEN snow(ki) = snow(ki)+precip_snow(ki)*dtime*(1.-snow_wfact*(1.-fsic(ki))) END IF ! snow and ice sublimation IF (evap(ki) > 0.) THEN snow_evap = MIN (snow(ki) / dtime, evap(ki)) snow(ki) = snow(ki) - snow_evap * dtime snow(ki) = MAX(0.0, snow(ki)) seaice(ki) = MAX(0.0,seaice(ki)-(evap(ki)-snow_evap)*dtime) ENDIF ! Melt / Freeze snow from above if Tsurf>0 IF (tsurf_new(ki).GT.t_melt) THEN ! energy available for melting snow (in kg of melted snow /m2) e_melt = MIN(MAX(snow(ki)*(tsurf_new(ki)-t_melt)*ice_cap/2. & /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(ki)) ! remove snow IF (snow(ki).GT.e_melt) THEN snow(ki)=snow(ki)-e_melt tsurf_new(ki)=t_melt ELSE ! all snow is melted ! add remaining heat flux to ice e_melt=e_melt-snow(ki) snow(ki)=0.0 tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki)) tsurf_new(ki)=tice(ki) END IF END IF ! Bottom melt / grow ! bottom freeze if bottom flux (cond + oce-ice) <0 IF (f_bot(ki).LT.0) THEN ! larger fraction of bottom growth frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thick) & / (h_ice_max-h_ice_thick))) ! quantity of new ice (formed at mean ice temp) e_melt= -f_bot(ki) * dtime * fsic(ki) & / (ice_lat+ice_cap/2.*(t_freeze-tice(ki))) ! first increase height to h_thick dhsic=MAX(0.,MIN(h_ice_thick-seaice(ki),e_melt/fsic(ki))) seaice(ki)=dhsic+seaice(ki) e_melt=e_melt-fsic(ki)*dhsic IF (e_melt.GT.0.) THEN ! frac_mf fraction used for lateral increase dfsic=MIN(ice_frac_max-fsic(ki),e_melt*frac_mf/seaice(ki)) fsic(ki)=fsic(ki)+dfsic e_melt=e_melt-dfsic*seaice(ki) ! rest used to increase height seaice(ki)=MIN(h_ice_max,seaice(ki)+e_melt/fsic(ki)) END IF ELSE ! melt from below if bottom flux >0 ! larger fraction of lateral melt from warm ocean frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thin) & / (h_ice_thick-h_ice_thin))) ! bring ice to freezing and melt from below ! quantity of melted ice e_melt= f_bot(ki) * dtime * fsic(ki) & / (ice_lat+ice_cap/2.*(tice(ki)-t_freeze)) ! first decrease height to h_thick IF (fsic(ki).GT.0) THEN dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki))) ELSE dhsic=0 ENDIF seaice(ki)=seaice(ki)-dhsic e_melt=e_melt-fsic(ki)*dhsic IF (e_melt.GT.0) THEN ! frac_mf fraction used for height decrease dhsic=MAX(0.,MIN(seaice(ki)-h_ice_min,e_melt*frac_mf/fsic(ki))) seaice(ki)=seaice(ki)-dhsic e_melt=e_melt-fsic(ki)*dhsic ! rest used to decrease fraction (up to 0!) dfsic=MIN(fsic(ki),e_melt/seaice(ki)) ! keep remaining in ocean if everything melted e_melt=e_melt-dfsic*seaice(ki) slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime ELSE dfsic=0 END IF fsic(ki)=fsic(ki)-dfsic END IF ! melt ice from above if Tice>0 IF (tice(ki).GT.t_melt) THEN ! quantity of ice melted (kg/m2) e_melt=MAX(seaice(ki)*(tice(ki)-t_melt)*ice_cap/2. & /(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0) ! melt from above, height only dhsic=MIN(seaice(ki)-h_ice_min,e_melt) e_melt=e_melt-dhsic IF (e_melt.GT.0) THEN ! lateral melt if ice too thin dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/h_ice_min*fsic(ki)) ! if all melted add remaining heat to ocean e_melt=MAX(0.,e_melt*fsic(ki)-dfsic*h_ice_min) slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime ! update height and fraction fsic(ki)=fsic(ki)-dfsic END IF seaice(ki)=seaice(ki)-dhsic ! surface temperature at melting point tice(ki)=t_melt tsurf_new(ki)=t_melt END IF ! convert snow to ice if below floating line h_test=(seaice(ki)+snow(ki))*ice_den-seaice(ki)*sea_den IF ((h_test.GT.0.).AND.(seaice(ki).GT.h_ice_min)) THEN !snow under water ! extra snow converted to ice (with added frozen sea water) dhsic=h_test/(sea_den-ice_den+sno_den) seaice(ki)=seaice(ki)+dhsic snow(ki)=snow(ki)-dhsic*sno_den/ice_den ! available energy (freeze sea water + bring to tice) e_melt=dhsic*(1.-sno_den/ice_den)*(ice_lat+ & ice_cap/2.*(t_freeze-tice(ki))) ! update ice temperature tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(ki)+seaice(ki)) END IF END DO !******************************************************************************* ! 3) cumulate ice-ocean fluxes, update tslab, lateral grow !***********************************************o******************************* !cumul fluxes. bilg_cum(:)=bilg_cum(:)+slab_bilg(:) IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab ! Add cumulated surface fluxes tslab(:,1)=tslab(:,1)+bilg_cum(:)*cyang*dtime bilg_cum(:)=0. ! If slab temperature below freezing, new lateral growth DO i=1,knon ki=knindex(i) IF (tslab(ki,1).LT.t_freeze) THEN ! create more ice ! quantity of new ice formed over open ocean ! (formed at mean ice temperature) e_freeze=(t_freeze-tslab(ki,1))/cyang & /(ice_lat+ice_cap/2.*(t_freeze-tice(ki))) ! new ice height and fraction h_new=MAX(h_ice_thin,MIN(h_ice_new,seaice(ki))) ! new height ! h_new=MIN(h_ice_new,seaice(ki)) dfsic=MIN(ice_frac_max-fsic(ki) & ,MAX(ice_frac_min,e_freeze/h_new)) ! update average sea ice height seaice(ki)=(seaice(ki)*fsic(ki)+e_freeze) & / (dfsic+fsic(ki)) ! update snow thickness snow(ki) = snow(ki) * fsic(ki) / (dfsic+fsic(ki)) ! update tslab to freezing tslab(ki,1)=t_freeze ! update sea ice fraction fsic(ki)=fsic(ki)+dfsic END IF ! tslab below freezing END DO END IF ! coupling time !**************************************************************************************** ! 4) Compute sea-ice and snow albedo !**************************************************************************************** ! Removed all albedo stuff as it is computed through hydrol in the Generic model ! ------ End Albedo ---------- !tests remaining ice fraction (on ocean grid points) WHERE ((zmasq==0).and. & ((fsic.LT.ice_frac_min).OR.(seaice.LT.h_ice_min))) !! WHERE ((fsic.LT.ice_frac_min).OR.(seaice.LT.h_ice_min)) tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang tice=t_melt fsic=0. seaice=0. END WHERE END SUBROUTINE ocean_slab_ice ! !**************************************************************************************** ! SUBROUTINE ocean_slab_final !**************************************************************************************** ! Deallocate module variables !**************************************************************************************** IF (ALLOCATED(tslab)) DEALLOCATE(tslab) IF (ALLOCATED(fsic)) DEALLOCATE(fsic) IF (ALLOCATED(tice)) DEALLOCATE(tice) IF (ALLOCATED(seaice)) DEALLOCATE(seaice) IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg) IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum) IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum) IF (ALLOCATED(taux_cum)) DEALLOCATE(taux_cum) IF (ALLOCATED(tauy_cum)) DEALLOCATE(tauy_cum) IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman) IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff) IF (ALLOCATED(dt_gm)) DEALLOCATE(dt_gm) ! IF (ALLOCATED(dt_qflux)) DEALLOCATE(dt_qflux) ! IF (ALLOCATED(dt_qflux_sic)) DEALLOCATE(dt_qflux_sic) END SUBROUTINE ocean_slab_final ! !**************************************************************************************** ! SUBROUTINE ocean_slab_get_vars(ngrid, tslab_loc, tice_loc, seaice_loc, & flux_g_loc, dt_hdiff_loc,dt_ekman_loc,dt_gm_loc) ! "Get some variables from module ocean_slab_mod" INTEGER, INTENT(IN) :: ngrid REAL, DIMENSION(ngrid,nslay), INTENT(OUT) :: tslab_loc REAL, INTENT(OUT) :: tice_loc(ngrid) REAL, DIMENSION(ngrid), INTENT(OUT) :: seaice_loc REAL, DIMENSION(ngrid), INTENT(OUT) :: flux_g_loc REAL, DIMENSION(ngrid,nslay), INTENT(OUT) :: dt_hdiff_loc ! [in W/m2] REAL, DIMENSION(ngrid,nslay), INTENT(OUT) :: dt_ekman_loc ! [in W/m2] REAL, DIMENSION(ngrid,nslay), INTENT(OUT) :: dt_gm_loc ! [in W/m2] INTEGER :: i ! Set the output variables tslab_loc(:,:) = 0. tice_loc(:)=0. dt_hdiff_loc(:,:)=0. dt_ekman_loc(:,:)=0. dt_gm_loc(:,:)=0. tslab_loc(:,:) = tslab(:,:) tice_loc(:)=tice(:) seaice_loc(:) = seaice(:) flux_g_loc(:) = slab_bilg(:) !! dt_hdiff_loc(:,i) = dt_hdiff(:,i)*slabh(i)*1000.*4228. !Convert en W/m2 !! dt_ekman_loc(:,i) = dt_ekman(:,i)*slabh(i)*1000.*4228. IF (slab_hdiff) THEN DO i=1,nslay dt_hdiff_loc(:,i) = dt_hdiff(:,i) ENDDO ENDIF IF (slab_ekman) THEN DO i=1,nslay dt_ekman_loc(:,i) = dt_ekman(:,i) ENDDO ENDIF IF (slab_gm) THEN DO i=1,nslay dt_gm_loc(:,i) = dt_gm(:,i) ENDDO ENDIF END SUBROUTINE ocean_slab_get_vars ! !**************************************************************************************** ! END MODULE ocean_slab_mod