Ignore:
Timestamp:
Jun 16, 2014, 4:33:38 PM (10 years ago)
Author:
Ehouarn Millour
Message:

Preparatory stuff to fix and improve the slab ocean model.
FC

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r1907 r2057  
    55! "ocean=slab".
    66!
     7
     8  USE dimphy
     9  USE indice_sol_mod
     10
    711  IMPLICIT NONE
    812  PRIVATE
    9   PUBLIC :: ocean_slab_frac, ocean_slab_noice
     13  PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice!, ocean_slab_ice
     14
     15  INTEGER, PRIVATE, SAVE                           :: cpl_pas
     16  !$OMP THREADPRIVATE(cpl_pas)
     17  REAL, PRIVATE, SAVE                              :: cyang
     18  !$OMP THREADPRIVATE(cyang)
     19  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slabh
     20  !$OMP THREADPRIVATE(slabh)
     21  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: tslab
     22  !$OMP THREADPRIVATE(tslab)
     23  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE  :: pctsrf
     24  !$OMP THREADPRIVATE(pctsrf)
     25  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bils
     26  !$OMP THREADPRIVATE(slab_bils)
     27  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
     28  !$OMP THREADPRIVATE(bils_cum)
    1029
    1130CONTAINS
    … …  
    1332!****************************************************************************************
    1433!
    15   SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf, is_modified)
    16 
    17     USE dimphy
     34  SUBROUTINE ocean_slab_init(dtime, pctsrf_rst)
     35  !, seaice_rst etc
     36
     37    use IOIPSL
     38
     39    INCLUDE "iniprint.h"
     40    ! For ok_xxx vars (Ekman...)
     41    INCLUDE "clesphys.h"
     42
     43    ! Input variables
     44!****************************************************************************************
     45    REAL, INTENT(IN)                         :: dtime
     46! Variables read from restart file
     47    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
     48
     49! Local variables
     50!****************************************************************************************
     51    INTEGER                :: error
     52    CHARACTER (len = 80)   :: abort_message
     53    CHARACTER (len = 20)   :: modname = 'ocean_slab_intit'
     54
     55!****************************************************************************************
     56! Allocate surface fraction read from restart file
     57!****************************************************************************************
     58    ALLOCATE(pctsrf(klon,nbsrf), stat = error)
     59    IF (error /= 0) THEN
     60       abort_message='Pb allocation tmp_pctsrf_slab'
     61       CALL abort_gcm(modname,abort_message,1)
     62    ENDIF
     63    pctsrf(:,:) = pctsrf_rst(:,:)
     64
     65!****************************************************************************************
     66! Allocate local variables
     67!****************************************************************************************
     68    ALLOCATE(slab_bils(klon), stat = error)
     69    IF (error /= 0) THEN
     70       abort_message='Pb allocation slab_bils'
     71       CALL abort_gcm(modname,abort_message,1)
     72    ENDIF
     73    slab_bils(:) = 0.0   
     74    ALLOCATE(bils_cum(klon), stat = error)
     75    IF (error /= 0) THEN
     76       abort_message='Pb allocation slab_bils_cum'
     77       CALL abort_gcm(modname,abort_message,1)
     78    ENDIF
     79    bils_cum(:) = 0.0   
     80
     81! Layer thickness
     82    ALLOCATE(slabh(nslay), stat = error)
     83    IF (error /= 0) THEN
     84       abort_message='Pb allocation slabh'
     85       CALL abort_gcm(modname,abort_message,1)
     86    ENDIF
     87    slabh(1)=50.
     88! cyang = 1/heat capacity of top layer (rho.c.H)
     89    cyang=1/(slabh(1)*4.228e+06)
     90
     91! cpl_pas periode de couplage avec slab (update tslab, pctsrf)
     92! pour un calcul à chaque pas de temps, cpl_pas=1
     93    cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
     94    CALL getin('cpl_pas',cpl_pas)
     95    print *,'cpl_pas',cpl_pas
     96 END SUBROUTINE ocean_slab_init
     97!
     98!****************************************************************************************
     99!
     100  SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified)
     101
    18102    USE limit_read_mod
    19103    USE surface_data
    20     USE indice_sol_mod
    21104
    22105!    INCLUDE "clesphys.h"
    … …  
    27110    INTEGER, INTENT(IN)                        :: jour    ! jour a lire dans l'annee
    28111    REAL   , INTENT(IN)                        :: dtime   ! pas de temps de la physique (en s)
    29     REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf  ! sub-surface fraction
     112    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg  ! sub-surface fraction
    30113    LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is modified at this time step
    31114
    32115! Local variables
    33116!****************************************************************************************
    34     CHARACTER (len = 80)   :: abort_message
    35     CHARACTER (len = 20)   :: modname = 'ocean_slab_frac'
    36 
    37 
    38     IF (version_ocean == 'sicOBS') THEN   
    39        CALL limit_read_frac(itime, dtime, jour, pctsrf, is_modified)
     117
     118    IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN   
     119       CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified)
    40120    ELSE
    41        abort_message='Ocean slab model without forced sea-ice fractions has to be rewritten!!!'
    42        CALL abort_gcm(modname,abort_message,1)
    43 ! Here should sea-ice/ocean fraction either be calculated or returned if saved as a module varaiable
    44 ! (in the case the new fractions are calculated in ocean_slab_ice or ocean_slab_noice subroutines). 
     121       pctsrf_chg(:,:)=pctsrf(:,:)
     122       is_modified=.TRUE.
    45123    END IF
    46124
    … …  
    57135       radsol, snow, agesno, &
    58136       qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
    59        tsurf_new, dflux_s, dflux_l, lmt_bils)
     137       tsurf_new, dflux_s, dflux_l, qflux)
    60138   
    61     USE dimphy
    62139    USE calcul_fluxs_mod
    63     USE indice_sol_mod
     140    USE surface_data
    64141
    65142    INCLUDE "iniprint.h"
    … …  
    95172    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new
    96173    REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l     
    97     REAL, DIMENSION(klon), INTENT(OUT)   :: lmt_bils
     174    REAL, DIMENSION(klon), INTENT(OUT)   :: qflux
    98175
    99176! Local variables
    100177!****************************************************************************************
    101     INTEGER               :: i
     178    INTEGER               :: i,ki
    102179    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
    103     REAL, DIMENSION(klon) :: lmt_bils_oce, lmt_foce, diff_sst
     180    REAL, DIMENSION(klon) :: diff_sst, lmt_bils
    104181    REAL, DIMENSION(klon) :: u0, v0
    105182    REAL, DIMENSION(klon) :: u1_lay, v1_lay
    106     REAL                  :: calc_bils_oce, deltat
    107     REAL, PARAMETER       :: cyang=50.0 * 4.228e+06 ! capacite calorifique volumetrique de l'eau J/(m2 K)
    108183
    109184!****************************************************************************************
    … …  
    136211         flux_u1, flux_v1) 
    137212
    138 !****************************************************************************************
    139 ! 2) Get global variables lmt_bils and lmt_foce from file limit_slab.nc
    140 !
    141 !****************************************************************************************
    142     CALL limit_slab(itime, dtime, jour, lmt_bils, lmt_foce, diff_sst)  ! global pour un processus
    143 
    144     lmt_bils_oce(:) = 0.
    145     WHERE (lmt_foce > 0.)
    146        lmt_bils_oce = lmt_bils / lmt_foce ! global
    147     END WHERE
     213! Accumulate total fluxes locally
     214    slab_bils(:)=0.
     215    DO i=1,knon
     216        ki=knindex(i)
     217        slab_bils(ki)=(fluxlat(i)+fluxsens(i)+radsol(i))*pctsrf(ki,is_oce)/(1.-zmasq(ki))
     218        bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
     219! Also taux, tauy, saved vars...
     220    END DO
     221
     222!****************************************************************************************
     223! 2) Get global variables lmt_bils and diff_sst from file limit_slab.nc
     224!
     225!****************************************************************************************
     226    lmt_bils(:)=0.
     227    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst)  ! global pour un processus
     228    ! lmt_bils and diff_sst saved by limit_slab
     229    qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.
     230    ! qflux = total QFlux correction (in W/m2)
    148231
    149232!****************************************************************************************
    150233! 3) Recalculate new temperature
    151234!
    152 !****************************************************************************************
    153     DO i = 1, knon
    154        calc_bils_oce = radsol(i) + fluxsens(i) + fluxlat(i)
    155        deltat        = (calc_bils_oce - lmt_bils_oce(knindex(i)))*dtime/cyang +diff_sst(knindex(i))
    156        tsurf_new(i)  = tsurf_in(i) + deltat
    157     END DO
    158 
     235!***********************************************o*****************************************
     236    tsurf_new=tsurf_in
     237    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
     238        ! Compute transport
     239        ! Add read QFlux and SST tendency
     240        tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas
     241        ! Add cumulated surface fluxes
     242        tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
     243        ! Update surface temperature
     244        SELECT CASE(version_ocean)
     245        CASE('sicNO')
     246            DO i=1,knon
     247                ki=knindex(i)
     248                tsurf_new(i)=tslab(ki,1)
     249            END DO
     250        CASE('sicOBS') ! check for sea ice or tsurf below freezing
     251            DO i=1,knon
     252                ki=knindex(i)
     253                IF ((tslab(ki,1).LT.t_freeze).OR.(pctsrf(ki,is_sic).GT.epsfra)) THEN
     254                    tsurf_new(i)=t_freeze
     255                    tslab(ki,1)=t_freeze
     256                ELSE
     257                    tsurf_new(i)=tslab(ki,1)
     258                END IF
     259            END DO
     260        CASE('sicINT')
     261            DO i=1,knon
     262                ki=knindex(i)
     263                IF (pctsrf(ki,is_sic).LT.epsfra) THEN ! Free of ice
     264                    IF (tslab(ki,1).GT.t_freeze) THEN
     265                        tsurf_new(i)=tslab(ki,1)
     266                    ELSE
     267                        tsurf_new(i)=t_freeze
     268                        ! Call new ice routine
     269                        tslab(ki,1)=t_freeze
     270                    END IF
     271                ELSE ! ice present, tslab update completed in slab_ice
     272                    tsurf_new(i)=t_freeze
     273                END IF !ice free
     274            END DO
     275        END SELECT
     276        bils_cum(:)=0.0! clear cumulated fluxes
     277    END IF ! coupling time
    159278  END SUBROUTINE ocean_slab_noice
    160279!
    161280!****************************************************************************************
    162281!
     282!  SUBROUTINE ocean_slab_ice(   &
     283!       itime, dtime, jour, knon, knindex, &
     284!       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
     285!       AcoefH, AcoefQ, BcoefH, BcoefQ, &
     286!       AcoefU, AcoefV, BcoefU, BcoefV, &
     287!       ps, u1, v1, &
     288!       radsol, snow, qsurf, qsol, agesno, tsoil, &
     289!       alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
     290!       tsurf_new, dflux_s, dflux_l)
     291!
     292!****************************************************************************************
     293! 1) Flux calculation
     294!****************************************************************************************
     295! set beta, cal etc. depends snow / ice surf ?
     296! calcul_fluxs (sens, lat etc)
     297! calcul_flux_wind
     298
     299!****************************************************************************************
     300! 2) Update surface
     301!****************************************************************************************
     302! neige, fonte
     303! flux glace-ocean
     304! update temperature
     305! neige precip, evap
     306! Melt snow & ice from above
     307! New albedo
     308
     309!****************************************************************************************
     310! 3) Recalculate new ocean temperature
     311!    Melt / freeze from below
     312!***********************************************o*****************************************
     313
     314
     315!  END SUBROUTINE ocean_slab_ice
     316!
     317!****************************************************************************************
     318!
     319  SUBROUTINE ocean_slab_final
     320  !, seaice_rst etc
     321
     322    ! For ok_xxx vars (Ekman...)
     323    INCLUDE "clesphys.h"
     324
     325!****************************************************************************************
     326! Deallocate module variables
     327!
     328!****************************************************************************************
     329    IF (ALLOCATED(pctsrf)) DEALLOCATE(pctsrf)
     330    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
     331
     332  END SUBROUTINE ocean_slab_final
     333
    163334END MODULE ocean_slab_mod
Note: See TracChangeset for help on using the changeset viewer.