Ignore:
Timestamp:
Oct 10, 2016, 10:57:24 AM (8 years ago)
Author:
Ehouarn Millour
Message:

Making the slab work:

  • added a slab_heat_transp_mod module for horizontal diffusion and Ekman transport
  • added storage and output of relevent variables in phyredem, phyetat0, phy_output_ctrlout_mod, phys_output_write_mod
  • moved nslay (number of slab layers) out of dimphy into ocean_slab_mod.

FC

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/ocean_slab_mod.F90

    r2311 r2656  
    99  USE indice_sol_mod
    1010  USE surface_data
     11  USE mod_grid_phy_lmdz, ONLY: klon_glo
     12  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
    1113
    1214  IMPLICIT NONE
    … …  
    1416  PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice
    1517
    16 !****************************************************************************************
     18!***********************************************************************************
    1719! Global saved variables
    18 !****************************************************************************************
    19 
     20!***********************************************************************************
     21  ! number of slab vertical layers
     22  INTEGER, PUBLIC, SAVE :: nslay
     23  !$OMP THREADPRIVATE(nslay)
     24  ! timestep for coupling (update slab temperature) in timesteps
    2025  INTEGER, PRIVATE, SAVE                           :: cpl_pas
    2126  !$OMP THREADPRIVATE(cpl_pas)
     27  ! cyang = 1/heat capacity of top layer (rho.c.H)
    2228  REAL, PRIVATE, SAVE                              :: cyang
    2329  !$OMP THREADPRIVATE(cyang)
     30  ! depth of slab layers (1 or 2)
    2431  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slabh
    2532  !$OMP THREADPRIVATE(slabh)
     33  ! slab temperature
    2634  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: tslab
    2735  !$OMP THREADPRIVATE(tslab)
     36  ! heat flux convergence due to Ekman
     37  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_ekman
     38  !$OMP THREADPRIVATE(dt_ekman)
     39  ! heat flux convergence due to horiz diffusion
     40  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_hdiff
     41  !$OMP THREADPRIVATE(dt_hdiff)
     42  ! fraction of ocean covered by sea ice (sic / (oce+sic))
    2843  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: fsic
    2944  !$OMP THREADPRIVATE(fsic)
     45  ! temperature of the sea ice
    3046  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: tice
    3147  !$OMP THREADPRIVATE(tice)
     48  ! sea ice thickness, in kg/m2
    3249  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: seaice
    3350  !$OMP THREADPRIVATE(seaice)
     51  ! net surface heat flux, weighted by open ocean fraction
    3452  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bils
    3553  !$OMP THREADPRIVATE(slab_bils)
     54  ! slab_bils accumulated over cpl_pas timesteps
    3655  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
    3756  !$OMP THREADPRIVATE(bils_cum)
     57  ! net heat flux into the ocean below the ice : conduction + solar radiation
    3858  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bilg
    3959  !$OMP THREADPRIVATE(slab_bilg)
     60  ! slab_bilg over cpl_pas timesteps
    4061  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bilg_cum
    4162  !$OMP THREADPRIVATE(bilg_cum)
    42 
    43 !****************************************************************************************
     63  ! wind stress saved over cpl_pas timesteps
     64  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: taux_cum
     65  !$OMP THREADPRIVATE(taux_cum)
     66  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: tauy_cum
     67  !$OMP THREADPRIVATE(tauy_cum)
     68
     69!***********************************************************************************
    4470! Parameters (could be read in def file: move to slab_init)
    45 !****************************************************************************************
     71!***********************************************************************************
    4672! snow and ice physical characteristics:
    4773    REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp
    4874    REAL, PARAMETER :: t_melt=273.15   ! melting ice temp
    4975    REAL, PARAMETER :: sno_den=300. !mean snow density, kg/m3
    50     REAL, PARAMETER :: ice_den=917.
    51     REAL, PARAMETER :: sea_den=1025.
    52     REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity
    53     REAL, PARAMETER :: sno_cond=0.31*sno_den
     76    REAL, PARAMETER :: ice_den=917. ! ice density
     77    REAL, PARAMETER :: sea_den=1025. ! sea water density
     78    REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity of ice
     79    REAL, PARAMETER :: sno_cond=0.31*sno_den ! conductivity of snow
    5480    REAL, PARAMETER :: ice_cap=2067.   ! specific heat capacity, snow and ice
     81    REAL, PARAMETER :: sea_cap=3995.   ! specific heat capacity, snow and ice
    5582    REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice
    5683
    … …  
    74101    REAL, PARAMETER :: alb_ice_dry=0.75 !dry thick ice
    75102    REAL, PARAMETER :: alb_ice_wet=0.66 !melting thick ice
    76     REAL, PARAMETER :: pen_frac=0.3 !fraction of penetrating shortwave (no snow)
     103    REAL, PARAMETER :: pen_frac=0.3 !fraction of shortwave penetrating into the
     104    ! ice (no snow)
    77105    REAL, PARAMETER :: pen_ext=1.5 !extinction of penetrating shortwave (m-1)
    78106
    79 !****************************************************************************************
     107! horizontal transport
     108   LOGICAL, PUBLIC, SAVE :: slab_hdiff
     109   !$OMP THREADPRIVATE(slab_hdiff)
     110   REAL, PRIVATE, SAVE    :: coef_hdiff ! coefficient for horizontal diffusion
     111   !$OMP THREADPRIVATE(coef_hdiff)
     112   INTEGER, PUBLIC, SAVE :: slab_ekman, slab_cadj ! Ekman, conv adjustment
     113   !$OMP THREADPRIVATE(slab_ekman)
     114   !$OMP THREADPRIVATE(slab_cadj)
     115
     116!***********************************************************************************
    80117
    81118CONTAINS
    82119!
    83 !****************************************************************************************
     120!***********************************************************************************
    84121!
    85122  SUBROUTINE ocean_slab_init(dtime, pctsrf_rst)
    86123  !, seaice_rst etc
    87124
    88     use IOIPSL
    89 
    90     ! For ok_xxx vars (Ekman...)
    91     INCLUDE "clesphys.h"
     125    USE ioipsl_getin_p_mod, ONLY : getin_p
     126    USE mod_phys_lmdz_transfert_para, ONLY : gather
     127    USE slab_heat_transp_mod, ONLY : ini_slab_transp
    92128
    93129    ! Input variables
    94 !****************************************************************************************
     130!***********************************************************************************
    95131    REAL, INTENT(IN)                         :: dtime
    96132! Variables read from restart file
    97     REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
     133    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
     134    ! surface fractions from start file
    98135
    99136! Local variables
    100 !****************************************************************************************
     137!************************************************************************************
    101138    INTEGER                :: error
     139    REAL, DIMENSION(klon_glo) :: zmasq_glo
    102140    CHARACTER (len = 80)   :: abort_message
    103141    CHARACTER (len = 20)   :: modname = 'ocean_slab_intit'
    104142
    105 !****************************************************************************************
     143!***********************************************************************************
     144! Define some parameters
     145!***********************************************************************************
     146! Number of slab layers
     147    nslay=2
     148    CALL getin_p('slab_layers',nslay)
     149    print *,'number of slab layers : ',nslay
     150! Layer thickness
     151    ALLOCATE(slabh(nslay), stat = error)
     152    IF (error /= 0) THEN
     153       abort_message='Pb allocation slabh'
     154       CALL abort_physic(modname,abort_message,1)
     155    ENDIF
     156    slabh(1)=50.
     157    IF (nslay.GT.1) THEN
     158        slabh(2)=150.
     159    END IF
     160
     161! cyang = 1/heat capacity of top layer (rho.c.H)
     162    cyang=1/(slabh(1)*sea_den*sea_cap)
     163
     164! cpl_pas  coupling period (update of tslab and ice fraction)
     165! pour un calcul a chaque pas de temps, cpl_pas=1
     166    cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
     167    CALL getin_p('cpl_pas',cpl_pas)
     168    print *,'cpl_pas',cpl_pas
     169
     170! Horizontal diffusion
     171    slab_hdiff=.FALSE.
     172    CALL getin_p('slab_hdiff',slab_hdiff)
     173    coef_hdiff=25000.
     174    CALL getin_p('coef_hdiff',coef_hdiff)
     175! Ekman transport
     176    slab_ekman=0
     177    CALL getin_p('slab_ekman',slab_ekman)
     178! Convective adjustment
     179    IF (nslay.EQ.1) THEN
     180        slab_cadj=0
     181    ELSE
     182        slab_cadj=1
     183    END IF
     184    CALL getin_p('slab_cadj',slab_cadj)
     185
     186!************************************************************************************
    106187! Allocate surface fraction read from restart file
    107 !****************************************************************************************
     188!************************************************************************************
    108189    ALLOCATE(fsic(klon), stat = error)
    109190    IF (error /= 0) THEN
    … …  
    112193    ENDIF
    113194    fsic(:)=0.
     195    !zmasq = continent fraction
    114196    WHERE (1.-zmasq(:)>EPSFRA)
    115197        fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:))
    116198    END WHERE
    117199
    118 !****************************************************************************************
    119 ! Allocate saved variables
    120 !****************************************************************************************
     200!************************************************************************************
     201! Allocate saved fields
     202!************************************************************************************
    121203    ALLOCATE(tslab(klon,nslay), stat=error)
    122204       IF (error /= 0) CALL abort_physic &
    … …  
    136218    bils_cum(:) = 0.0   
    137219
    138     IF (version_ocean=='sicINT') THEN
     220    IF (version_ocean=='sicINT') THEN ! interactive sea ice
    139221        ALLOCATE(slab_bilg(klon), stat = error)
    140222        IF (error /= 0) THEN
    … …  
    161243    END IF
    162244
    163 !****************************************************************************************
    164 ! Define some parameters
    165 !****************************************************************************************
    166 ! Layer thickness
    167     ALLOCATE(slabh(nslay), stat = error)
    168     IF (error /= 0) THEN
    169        abort_message='Pb allocation slabh'
    170        CALL abort_physic(modname,abort_message,1)
     245    IF (slab_hdiff) THEN !horizontal diffusion
     246        ALLOCATE(dt_hdiff(klon,nslay), stat = error)
     247        IF (error /= 0) THEN
     248           abort_message='Pb allocation dt_hdiff'
     249           CALL abort_physic(modname,abort_message,1)
     250        ENDIF
     251        dt_hdiff(:,:) = 0.0   
    171252    ENDIF
    172     slabh(1)=50.
    173 ! cyang = 1/heat capacity of top layer (rho.c.H)
    174     cyang=1/(slabh(1)*4.228e+06)
    175 
    176 ! cpl_pas periode de couplage avec slab (update tslab, pctsrf)
    177 ! pour un calcul à chaque pas de temps, cpl_pas=1
    178     cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
    179     CALL getin('cpl_pas',cpl_pas)
    180     print *,'cpl_pas',cpl_pas
     253
     254    IF (slab_ekman.GT.0) THEN ! ekman transport
     255        ALLOCATE(dt_ekman(klon,nslay), stat = error)
     256        IF (error /= 0) THEN
     257           abort_message='Pb allocation dt_ekman'
     258           CALL abort_physic(modname,abort_message,1)
     259        ENDIF
     260        dt_ekman(:,:) = 0.0   
     261        ALLOCATE(taux_cum(klon), stat = error)
     262        IF (error /= 0) THEN
     263           abort_message='Pb allocation taux_cum'
     264           CALL abort_physic(modname,abort_message,1)
     265        ENDIF
     266        taux_cum(:) = 0.0   
     267        ALLOCATE(tauy_cum(klon), stat = error)
     268        IF (error /= 0) THEN
     269           abort_message='Pb allocation tauy_cum'
     270           CALL abort_physic(modname,abort_message,1)
     271        ENDIF
     272        tauy_cum(:) = 0.0   
     273    ENDIF
     274
     275! Initialize transport
     276    IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN
     277      CALL gather(zmasq,zmasq_glo)
     278!$OMP MASTER  ! Only master thread
     279      IF (is_mpi_root) THEN
     280          CALL ini_slab_transp(zmasq_glo)
     281      END IF
     282!$OMP END MASTER
     283    END IF
    181284
    182285 END SUBROUTINE ocean_slab_init
    183286!
    184 !****************************************************************************************
     287!***********************************************************************************
    185288!
    186289  SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified)
    187290
    188     USE limit_read_mod
    189 
    190 !    INCLUDE "clesphys.h"
     291! this routine sends back the sea ice and ocean fraction to the main physics
     292! routine. Called only with interactive sea ice
    191293
    192294! Arguments
    193 !****************************************************************************************
    194     INTEGER, INTENT(IN)                        :: itime   ! numero du pas de temps courant
    195     INTEGER, INTENT(IN)                        :: jour    ! jour a lire dans l'annee
    196     REAL   , INTENT(IN)                        :: dtime   ! pas de temps de la physique (en s)
     295!************************************************************************************
     296    INTEGER, INTENT(IN)                        :: itime   ! current timestep
     297    INTEGER, INTENT(IN)                        :: jour    !  day in year (not
     298    REAL   , INTENT(IN)                        :: dtime   ! physics timestep (s)
    197299    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg  ! sub-surface fraction
    198     LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is modified at this time step
    199 
    200 ! Local variables
    201 !****************************************************************************************
    202 
    203     IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN   
    204        CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified)
    205     ELSE
     300    LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is
     301                                                         ! modified at this time step
     302
    206303       pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
    207304       pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:))
    208305       is_modified=.TRUE.
    209     END IF
    210306
    211307  END SUBROUTINE ocean_slab_frac
    212308!
    213 !****************************************************************************************
     309!************************************************************************************
    214310!
    215311  SUBROUTINE ocean_slab_noice( &
    … …  
    224320   
    225321    USE calcul_fluxs_mod
     322    USE slab_heat_transp_mod, ONLY: divgrad_phy, slab_ekman1, slab_ekman2
     323    USE mod_phys_lmdz_para
    226324
    227325    INCLUDE "clesphys.h"
    228326
     327! This routine
     328! (1) computes surface turbulent fluxes over points with some open ocean   
     329! (2) reads additional Q-flux (everywhere)
     330! (3) computes horizontal transport (diffusion & Ekman)
     331! (4) updates slab temperature every cpl_pas ; creates new ice if needed.
     332
     333! Note :
     334! klon total number of points
     335! knon number of points with open ocean (varies with time)
     336! knindex gives position of the knon points within klon.
     337! In general, local saved variables have klon values
     338! variables exchanged with PBL module have knon.
     339
    229340! Input arguments
    230 !****************************************************************************************
    231     INTEGER, INTENT(IN)                  :: itime
    232     INTEGER, INTENT(IN)                  :: jour
    233     INTEGER, INTENT(IN)                  :: knon
    234     INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
    235     REAL, INTENT(IN)                     :: dtime
    236     REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
    237     REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragq, cdragm
    238     REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
    239     REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
    240     REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
    241     REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
    242     REAL, DIMENSION(klon), INTENT(IN)    :: ps
    243     REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness
    244     REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
    245     REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
     341!***********************************************************************************
     342    INTEGER, INTENT(IN)                  :: itime ! current timestep INTEGER,
     343    INTEGER, INTENT(IN)                  :: jour  ! day in year (for Q-Flux)
     344    INTEGER, INTENT(IN)                  :: knon  ! number of points
     345    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
     346    REAL, INTENT(IN) :: dtime  ! timestep (s)
     347    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
     348    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragq, cdragm
     349    ! drag coefficients
     350    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
     351    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum ! near surface T, q
     352    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
     353    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
     354    ! exchange coefficients for boundary layer scheme
     355    REAL, DIMENSION(klon), INTENT(IN)    :: ps  ! surface pressure
     356    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness ! surface wind
     357    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in ! surface temperature
     358    REAL, DIMENSION(klon), INTENT(INOUT) :: radsol ! net surface radiative flux
    246359
    247360! In/Output arguments
    248 !****************************************************************************************
    249     REAL, DIMENSION(klon), INTENT(INOUT) :: snow
     361!************************************************************************************
     362    REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! in kg/m2
    250363   
    251364! Output arguments
    252 !****************************************************************************************
     365!************************************************************************************
    253366    REAL, DIMENSION(klon), INTENT(OUT)   :: qsurf
    254367    REAL, DIMENSION(klon), INTENT(OUT)   :: evap, fluxsens, fluxlat
    255368    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u1, flux_v1
    256     REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new
     369    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new ! new surface tempearture
    257370    REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l     
    258371    REAL, DIMENSION(klon), INTENT(OUT)   :: qflux
    259372
    260373! Local variables
    261 !****************************************************************************************
    262     INTEGER               :: i,ki
    263     ! surface fluxes
     374!************************************************************************************
     375    INTEGER               :: i,ki,k
     376    REAL                  :: t_cadj
     377    !  for surface heat fluxes
    264378    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
    265     ! flux correction
     379    ! for Q-Flux computation: d/dt SST, d/dt ice volume (kg/m2), surf fluxes
    266380    REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
    267     ! surface wind stress
     381    ! for surface wind stress
    268382    REAL, DIMENSION(klon) :: u0, v0
    269383    REAL, DIMENSION(klon) :: u1_lay, v1_lay
    270     ! ice creation
     384    ! for new ice creation
    271385    REAL                  :: e_freeze, h_new, dfsic
    272 
    273 !****************************************************************************************
    274 ! 1) Flux calculation
    275 !
    276 !****************************************************************************************
    277     cal(:)      = 0.
    278     beta(:)     = 1.
    279     dif_grnd(:) = 0.
     386    ! horizontal diffusion and Ekman local vars
     387    ! dimension = global domain (klon_glo) instead of // subdomains
     388    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo, dt_ekman_glo
     389    ! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop
     390    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_tmp, dt_ekman_tmp
     391    REAL, DIMENSION(klon_glo,nslay) :: tslab_glo
     392    REAL, DIMENSION(klon_glo) :: taux_glo,tauy_glo
     393
     394!****************************************************************************************
     395! 1) Surface fluxes calculation
     396!   
     397!****************************************************************************************
     398    cal(:)      = 0. ! infinite thermal inertia
     399    beta(:)     = 1. ! wet surface
     400    dif_grnd(:) = 0. ! no diffusion into ground
    280401   
    281402! Suppose zero surface speed
    … …  
    285406    v1_lay(:) = v1(:) - v0(:)
    286407
     408! Compute latent & sensible fluxes
    287409    CALL calcul_fluxs(knon, is_oce, dtime, &
    288410         tsurf_in, p1lay, cal, beta, cdragh, cdragq, ps, &
    … …  
    292414         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
    293415
    294 ! - Flux calculation at first modele level for U and V
    295     CALL calcul_flux_wind(knon, dtime, &
    296          u0, v0, u1, v1, gustiness, cdragm, &
    297          AcoefU, AcoefV, BcoefU, BcoefV, &
    298          p1lay, temp_air, &
    299          flux_u1, flux_v1) 
    300 
    301 ! Accumulate total fluxes locally
     416! save total cumulated heat fluxes locally
     417! radiative + turbulent + melt of falling snow
    302418    slab_bils(:)=0.
    303419    DO i=1,knon
    … …  
    306422                      -precip_snow(i)*ice_lat*(1.+snow_wfact*fsic(ki)))
    307423        bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
    308 ! Also taux, tauy, saved vars...
    309424    END DO
    310425
    311 !****************************************************************************************
    312 ! 2) Get global variables lmt_bils and diff_sst from file limit_slab.nc
     426!  Compute surface wind stress
     427    CALL calcul_flux_wind(knon, dtime, &
     428         u0, v0, u1, v1, gustiness, cdragm, &
     429         AcoefU, AcoefV, BcoefU, BcoefV, &
     430         p1lay, temp_air, &
     431         flux_u1, flux_v1) 
     432
     433! save cumulated wind stress
     434    IF (slab_ekman.GT.0) THEN
     435      DO i=1,knon
     436          ki=knindex(i)
     437          taux_cum(ki)=taux_cum(ki)+flux_u1(i)*(1.-fsic(ki))/cpl_pas
     438          tauy_cum(ki)=tauy_cum(ki)+flux_v1(i)*(1.-fsic(ki))/cpl_pas
     439      END DO
     440    ENDIF
     441
     442!****************************************************************************************
     443! 2) Q-Flux : get global variables lmt_bils, diff_sst and diff_siv from file limit_slab.nc
    313444!
    314445!****************************************************************************************
    315446    lmt_bils(:)=0.
    316     CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) ! global pour un processus
     447    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv)
    317448    ! lmt_bils and diff_sst,siv saved by limit_slab
    318449    qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400.
    … …  
    320451
    321452!****************************************************************************************
    322 ! 3) Recalculate new temperature
    323 !
     453! 3) Recalculate new temperature (add Surf fluxes, Q-Flux, Ocean transport)
     454!    Bring to freezing temp and make sea ice if necessary
     455! 
    324456!***********************************************o*****************************************
    325457    tsurf_new=tsurf_in
    326458    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
    327         ! Compute transport
    328         ! Add read QFlux and SST tendency
    329         tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas
    330         ! Add cumulated surface fluxes
    331         tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
    332         ! Update surface temperature
    333         SELECT CASE(version_ocean)
    334         CASE('sicNO')
    335             DO i=1,knon
    336                 ki=knindex(i)
    337                 tsurf_new(i)=tslab(ki,1)
     459! ***********************************
     460!  Horizontal transport
     461! ***********************************
     462      IF (slab_ekman.GT.0) THEN
     463          ! copy wind stress to global var
     464          CALL gather(taux_cum,taux_glo)
     465          CALL gather(tauy_cum,tauy_glo)
     466      END IF
     467
     468      IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN
     469        CALL gather(tslab,tslab_glo)
     470      ! Compute horiz transport on one process only
     471!$OMP MASTER  ! Only master thread
     472        IF (is_mpi_root) THEN ! Only master processus         
     473          IF (slab_hdiff) THEN
     474              dt_hdiff_glo(:,:)=0.
     475          END IF
     476          IF (slab_ekman.GT.0) THEN
     477              dt_ekman_glo(:,:)=0.
     478          END IF
     479          DO i=1,cpl_pas ! time splitting for numerical stability
     480            IF (slab_ekman.GT.0) THEN
     481              SELECT CASE (slab_ekman)
     482                CASE (1)
     483                  CALL slab_ekman1(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
     484                CASE (2)
     485                  CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
     486                CASE DEFAULT
     487                  dt_ekman_tmp(:,:)=0.
     488              END SELECT
     489              dt_ekman_glo(:,:)=dt_ekman_glo(:,:)+dt_ekman_tmp(:,:)
     490              ! convert dt_ekman from K.s-1.(kg.m-2) to K.s-1   
     491              DO k=1,nslay
     492                dt_ekman_tmp(:,k)=dt_ekman_tmp(:,k)/(slabh(k)*sea_den)
     493              ENDDO
     494              tslab_glo=tslab_glo+dt_ekman_tmp*dtime
     495            ENDIF
     496            IF (slab_hdiff) THEN ! horizontal diffusion
     497              ! laplacian of slab T
     498              CALL divgrad_phy(nslay,tslab_glo,dt_hdiff_tmp)
     499              ! multiply by diff coef and normalize to 50m slab equivalent
     500              dt_hdiff_tmp=dt_hdiff_tmp*coef_hdiff*50./SUM(slabh)
     501              dt_hdiff_glo(:,:)=dt_hdiff_glo(:,:)+ dt_hdiff_tmp(:,:)
     502              tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
     503            END IF
     504          END DO ! time splitting
     505          IF (slab_hdiff) THEN
     506            !dt_hdiff_glo saved in W/m2
     507            DO k=1,nslay
     508              dt_hdiff_glo(:,k)=dt_hdiff_glo(:,k)*slabh(k)*sea_den*sea_cap/cpl_pas
    338509            END DO
    339         CASE('sicOBS') ! check for sea ice or tslab below freezing
    340             DO i=1,knon
    341                 ki=knindex(i)
    342                 IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
    343                     tslab(ki,1)=t_freeze
    344                 END IF
    345                 tsurf_new(i)=tslab(ki,1)
    346             END DO
    347         CASE('sicINT')
    348             DO i=1,knon
    349                 ki=knindex(i)
    350                 IF (fsic(ki).LT.epsfra) THEN ! Free of ice
    351                     IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
    352                         ! quantity of new ice formed
    353                         e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
    354                         ! new ice
    355                         tice(ki)=t_freeze
    356                         fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
    357                         IF (fsic(ki).GT.ice_frac_min) THEN
    358                             seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
    359                             tslab(ki,1)=t_freeze
    360                         ELSE
    361                             fsic(ki)=0.
    362                         END IF
    363                         tsurf_new(i)=t_freeze
    364                     ELSE
    365                         tsurf_new(i)=tslab(ki,1)
    366                     END IF
    367                 ELSE ! ice present
    368                     tsurf_new(i)=t_freeze
    369                     IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
    370                         ! quantity of new ice formed over open ocean
    371                         e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
    372                                  /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
    373                         ! new ice height and fraction
    374                         h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
    375                         dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
    376                         h_new=MIN(e_freeze/dfsic,h_ice_max)
    377                         ! update tslab to freezing over open ocean only
    378                         tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki))
    379                         ! update sea ice
    380                         seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) &
    381                                    /(dfsic+fsic(ki))
    382                         fsic(ki)=fsic(ki)+dfsic
    383                         ! update snow?
    384                     END IF !freezing
    385                 END IF ! sea ice present
    386             END DO
    387         END SELECT
    388         bils_cum(:)=0.0! clear cumulated fluxes
     510          END IF
     511          IF (slab_ekman.GT.0) THEN
     512            ! dt_ekman_glo saved in W/m2
     513            dt_ekman_glo(:,:)=dt_ekman_glo(:,:)*sea_cap/cpl_pas
     514          END IF
     515        END IF ! is_mpi_root
     516!$OMP END MASTER
     517!$OMP BARRIER
     518        ! Send new fields back to all processes
     519        CALL Scatter(tslab_glo,tslab)
     520        IF (slab_hdiff) THEN
     521          CALL Scatter(dt_hdiff_glo,dt_hdiff)
     522        END IF
     523        IF (slab_ekman.GT.0) THEN
     524          CALL Scatter(dt_ekman_glo,dt_ekman)
     525          ! clear wind stress
     526          taux_cum(:)=0.
     527          tauy_cum(:)=0.
     528        END IF
     529      ENDIF ! transport
     530
     531! ***********************************
     532! Other heat fluxes
     533! ***********************************
     534      ! Add read QFlux
     535      tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas
     536      ! Add cumulated surface fluxes
     537      tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
     538      ! Convective adjustment if 2 layers
     539      IF ((nslay.GT.1).AND.(slab_cadj.GT.0)) THEN
     540        DO i=1,klon
     541          IF (tslab(i,2).GT.tslab(i,1)) THEN
     542            ! mean (mass-weighted) temperature
     543            t_cadj=SUM(tslab(i,:)*slabh(:))/SUM(slabh(:))
     544            tslab(i,1)=t_cadj
     545            tslab(i,2)=t_cadj
     546          END IF
     547        END DO
     548      END IF
     549! ***********************************
     550! Update surface temperature and ice
     551! ***********************************
     552      SELECT CASE(version_ocean)
     553      CASE('sicNO') ! no sea ice even below freezing !
     554          DO i=1,knon
     555              ki=knindex(i)
     556              tsurf_new(i)=tslab(ki,1)
     557          END DO
     558      CASE('sicOBS') ! "realistic" case, for prescribed sea ice
     559        ! tslab cannot be below freezing, or above it if there is sea ice
     560          DO i=1,knon
     561              ki=knindex(i)
     562              IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
     563                  tslab(ki,1)=t_freeze
     564              END IF
     565              tsurf_new(i)=tslab(ki,1)
     566          END DO
     567      CASE('sicINT') ! interactive sea ice
     568          DO i=1,knon
     569              ki=knindex(i)
     570              IF (fsic(ki).LT.epsfra) THEN ! Free of ice
     571                  IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
     572                      ! quantity of new ice formed
     573                      e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
     574                      ! new ice
     575                      tice(ki)=t_freeze
     576                      fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
     577                      IF (fsic(ki).GT.ice_frac_min) THEN
     578                          seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
     579                          tslab(ki,1)=t_freeze
     580                      ELSE
     581                          fsic(ki)=0.
     582                      END IF
     583                      tsurf_new(i)=t_freeze
     584                  ELSE
     585                      tsurf_new(i)=tslab(ki,1)
     586                  END IF
     587              ELSE ! ice present
     588                  tsurf_new(i)=t_freeze
     589                  IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
     590                      ! quantity of new ice formed over open ocean
     591                      e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
     592                               /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
     593                      ! new ice height and fraction
     594                      h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
     595                      dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
     596                      h_new=MIN(e_freeze/dfsic,h_ice_max)
     597                      ! update tslab to freezing over open ocean only
     598                      tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki))
     599                      ! update sea ice
     600                      seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) &
     601                                 /(dfsic+fsic(ki))
     602                      fsic(ki)=fsic(ki)+dfsic
     603                      ! update snow?
     604                  END IF ! tslab below freezing
     605              END IF ! sea ice present
     606          END DO
     607      END SELECT
     608      bils_cum(:)=0.0! clear cumulated fluxes
    389609    END IF ! coupling time
    390610  END SUBROUTINE ocean_slab_noice
    … …  
    530750           seaice(ki) = MAX(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime)
    531751        ENDIF
    532 ! Melt / Freeze from above if Tsurf>0
     752! Melt / Freeze snow from above if Tsurf>0
    533753        IF (tsurf_new(i).GT.t_melt) THEN
    534             ! energy available for melting snow (in kg/m2 of snow)
     754            ! energy available for melting snow (in kg of melted snow /m2)
    535755            e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. &
    536756               /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
    … …  
    672892    END IF ! coupling time
    673893   
    674     !tests fraction glace
     894    !tests ice fraction
    675895    WHERE (fsic.LT.ice_frac_min)
    676896        tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang
    … …  
    691911    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
    692912    IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
     913    IF (ALLOCATED(tice)) DEALLOCATE(tice)
     914    IF (ALLOCATED(seaice)) DEALLOCATE(seaice)
    693915    IF (ALLOCATED(slab_bils)) DEALLOCATE(slab_bils)
    694916    IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
    695917    IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
    696918    IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum)
    697     IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
     919    IF (ALLOCATED(taux_cum)) DEALLOCATE(taux_cum)
     920    IF (ALLOCATED(tauy_cum)) DEALLOCATE(tauy_cum)
     921    IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman)
     922    IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff)
    698923
    699924  END SUBROUTINE ocean_slab_final
Note: See TracChangeset for help on using the changeset viewer.