Changeset 3002


Ignore:
Timestamp:
Oct 3, 2017, 8:18:41 AM (7 years ago)
Author:
Ehouarn Millour
Message:

Improved slab routines.
FC

Location:
LMDZ5/trunk/libf
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90

    r2786 r3002  
    131131                                  cu,cuvsurcv,cv,cvusurcu, &
    132132                                  aire,apoln,apols, &
    133                                   aireu,airev)
     133                                  aireu,airev,rlatvdyn)
    134134  END IF
    135135

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

    r2989 r3002  
    24192419
    24202420    !--test on ocean surface albedo
    2421     IF (iflag_albedo.LT.0.OR.iflag_albedo.GT.1) THEN
     2421    IF (iflag_albedo.LT.0.OR.iflag_albedo.GT.2) THEN
    24222422       WRITE(lunout,*) ' ERROR iflag_albedo<>0,1'
    24232423       CALL abort_physic('conf_phys','choice iflag_albedo not valid',1)

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

    r2656 r3002  
    88  USE netcdf
    99  USE indice_sol_mod
     10  USE ocean_slab_mod, ONLY: nslay
    1011
    1112  IMPLICIT NONE
    … …  
    1819  INTEGER, INTENT(IN) :: jour    ! jour a lire dans l'annee
    1920  REAL   , INTENT(IN) :: dtime   ! pas de temps de la physique (en s)
    20   REAL, DIMENSION(klon), INTENT(OUT) :: lmt_bils, diff_sst, diff_siv
     21  REAL, DIMENSION(klon), INTENT(OUT) ::  diff_sst, diff_siv
     22  REAL, DIMENSION(klon,nslay), INTENT(OUT) :: lmt_bils
    2123
    2224! Locals variables with attribute SAVE
    2325!****************************************************************************************
    24   REAL, DIMENSION(:), ALLOCATABLE, SAVE :: bils_save, diff_sst_save
    25   REAL, DIMENSION(:), ALLOCATABLE, SAVE :: diff_siv_save
     26  REAL, DIMENSION(:), ALLOCATABLE, SAVE :: diff_siv_save, diff_sst_save
     27  REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: bils_save
    2628!$OMP THREADPRIVATE(bils_save, diff_sst_save, diff_siv_save)
    2729
    … …  
    3133  INTEGER                  :: nvarid, nid, ierr, i
    3234  INTEGER, DIMENSION(2)    :: start, epais
    33   REAL, DIMENSION(klon_glo):: bils_glo, sst_l_glo, sst_lp1_glo, diff_sst_glo
     35  REAL, DIMENSION(klon_glo):: sst_l_glo, sst_lp1_glo, diff_sst_glo
    3436  REAL, DIMENSION(klon_glo):: siv_l_glo, siv_lp1_glo, diff_siv_glo
     37  REAL, DIMENSION(klon_glo,nslay):: bils_glo
    3538  CHARACTER (len = 20)     :: modname = 'limit_slab'
     39  CHARACTER*2 str2
    3640  LOGICAL                  :: read_bils,read_sst,read_siv
    3741
    … …  
    4448! Initialize saved variables
    4549     IF (.NOT. ALLOCATED(bils_save)) THEN
    46         ALLOCATE(bils_save(klon), diff_sst_save(klon), diff_siv_save(klon), stat=ierr)
     50        ALLOCATE(bils_save(klon,nslay), diff_sst_save(klon), diff_siv_save(klon), stat=ierr)
    4751        IF (ierr /= 0) CALL abort_physic('limit_slab', 'pb in allocation',1)
    4852     END IF
    … …  
    7781!
    7882! Read bils_glo
    79         ierr = NF90_INQ_VARID(nid, 'BILS_OCE', nvarid)
     83        bils_glo(:,:)=0.
     84        ! First read first layer
     85        ! try first "BILS_OCE01"
     86        ierr = NF90_INQ_VARID(nid, 'BILS_OCE01', nvarid)
    8087        IF (ierr /= NF90_NOERR) THEN
    81             read_bils=.FALSE.
     88            ! Else BILS_OCE
     89            ierr = NF90_INQ_VARID(nid, 'BILS_OCE', nvarid)
     90            IF (ierr /= NF90_NOERR) THEN
     91              read_bils=.FALSE.
     92            ELSE
     93              ierr = NF90_GET_VAR(nid,nvarid,bils_glo(:,1),start,epais)
     94              IF (ierr /= NF90_NOERR) read_bils=.FALSE.
     95            ENDIF
    8296        ELSE
    83             ierr = NF90_GET_VAR(nid,nvarid,bils_glo,start,epais)
     97            ierr = NF90_GET_VAR(nid,nvarid,bils_glo(:,1),start,epais)
    8498            IF (ierr /= NF90_NOERR) read_bils=.FALSE.
    8599        END IF
     100        ! Try next layers if more than 1
     101        IF ((nslay.GT.1).AND.read_bils) THEN
     102          DO i=2,nslay
     103            WRITE(str2,'(i2.2)') i
     104            ierr = NF90_INQ_VARID(nid,'BILS_OCE'//str2, nvarid)
     105            IF (ierr.EQ.NF90_NOERR) THEN
     106              ierr = NF90_GET_VAR(nid,nvarid,bils_glo(:,i),start,epais)
     107            ENDIF
     108            IF (ierr /= NF90_NOERR) THEN
     109              print *,'WARNING : BILS_OCE not found for layer 2'
     110            ENDIF
     111          ENDDO
     112        ENDIF
     113
    86114! Read sst_glo for this day
    87115        ierr = NF90_INQ_VARID(nid, 'SST', nvarid)
    … …  
    146174         CALL Scatter(bils_glo, bils_save)
    147175     ELSE
    148          bils_save(:)=0.
     176         bils_save(:,:)=0.
    149177     END IF
    150178     IF (read_sst) THEN
    … …  
    161189  ENDIF ! time to read
    162190
    163   lmt_bils(:) = bils_save(:)
     191  lmt_bils(:,:) = bils_save(:,:)
    164192  diff_sst(:) = diff_sst_save(:)
    165193  diff_siv(:) = diff_siv_save(:)

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

    r2656 r3002  
    1 !
     1!Completed
    22MODULE ocean_slab_mod
    33!
    … …  
    4040  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_hdiff
    4141  !$OMP THREADPRIVATE(dt_hdiff)
     42  ! heat flux convergence due to GM eddy advection
     43  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_gm
     44  !$OMP THREADPRIVATE(dt_gm)
     45  ! Heat Flux correction
     46  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_qflux
     47  !$OMP THREADPRIVATE(dt_qflux)
    4248  ! fraction of ocean covered by sea ice (sic / (oce+sic))
    4349  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: fsic
    … …  
    5056  !$OMP THREADPRIVATE(seaice)
    5157  ! net surface heat flux, weighted by open ocean fraction
    52   REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bils
    53   !$OMP THREADPRIVATE(slab_bils)
    5458  ! slab_bils accumulated over cpl_pas timesteps
    5559  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
    … …  
    108112   LOGICAL, PUBLIC, SAVE :: slab_hdiff
    109113   !$OMP THREADPRIVATE(slab_hdiff)
     114   LOGICAL, PUBLIC, SAVE :: slab_gm
     115   !$OMP THREADPRIVATE(slab_gm)
    110116   REAL, PRIVATE, SAVE    :: coef_hdiff ! coefficient for horizontal diffusion
    111117   !$OMP THREADPRIVATE(coef_hdiff)
    … …  
    155161    ENDIF
    156162    slabh(1)=50.
     163    CALL getin_p('slab_depth',slabh(1))
    157164    IF (nslay.GT.1) THEN
    158165        slabh(2)=150.
    … …  
    176183    slab_ekman=0
    177184    CALL getin_p('slab_ekman',slab_ekman)
     185! GM eddy advection (2-layers only)
     186    slab_gm=.FALSE.
     187    CALL getin_p('slab_gm',slab_gm)
     188    IF (slab_ekman.LT.2) THEN
     189       slab_gm=.FALSE.
     190    ENDIF
    178191! Convective adjustment
    179192    IF (nslay.EQ.1) THEN
    … …  
    205218         (modname,'pb allocation tslab', 1)
    206219
    207     ALLOCATE(slab_bils(klon), stat = error)
    208     IF (error /= 0) THEN
    209        abort_message='Pb allocation slab_bils'
    210        CALL abort_physic(modname,abort_message,1)
    211     ENDIF
    212     slab_bils(:) = 0.0   
    213220    ALLOCATE(bils_cum(klon), stat = error)
    214221    IF (error /= 0) THEN
    … …  
    252259    ENDIF
    253260
     261    ALLOCATE(dt_qflux(klon,nslay), stat = error)
     262    IF (error /= 0) THEN
     263       abort_message='Pb allocation dt_qflux'
     264       CALL abort_physic(modname,abort_message,1)
     265    ENDIF
     266    dt_qflux(:,:) = 0.0   
     267
     268    IF (slab_gm) THEN !GM advection
     269        ALLOCATE(dt_gm(klon,nslay), stat = error)
     270        IF (error /= 0) THEN
     271           abort_message='Pb allocation dt_gm'
     272           CALL abort_physic(modname,abort_message,1)
     273        ENDIF
     274        dt_gm(:,:) = 0.0   
     275    ENDIF
     276
    254277    IF (slab_ekman.GT.0) THEN ! ekman transport
    255278        ALLOCATE(dt_ekman(klon,nslay), stat = error)
    … …  
    276299    IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN
    277300      CALL gather(zmasq,zmasq_glo)
    278 !$OMP MASTER  ! Only master thread
     301! Master thread/process only
     302!$OMP MASTER 
    279303      IF (is_mpi_root) THEN
    280304          CALL ini_slab_transp(zmasq_glo)
    … …  
    317341       radsol, snow, &
    318342       qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
    319        tsurf_new, dflux_s, dflux_l, qflux)
     343       tsurf_new, dflux_s, dflux_l, slab_bils)
    320344   
    321345    USE calcul_fluxs_mod
    322     USE slab_heat_transp_mod, ONLY: divgrad_phy, slab_ekman1, slab_ekman2
     346    USE slab_heat_transp_mod, ONLY: divgrad_phy,slab_ekman1,slab_ekman2,slab_gmdiff
    323347    USE mod_phys_lmdz_para
    324348
    … …  
    369393    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new ! new surface tempearture
    370394    REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l     
    371     REAL, DIMENSION(klon), INTENT(OUT)   :: qflux
     395    REAL, DIMENSION(klon), INTENT(OUT)   :: slab_bils
    372396
    373397! Local variables
    … …  
    378402    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
    379403    ! for Q-Flux computation: d/dt SST, d/dt ice volume (kg/m2), surf fluxes
    380     REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
     404    REAL, DIMENSION(klon) :: diff_sst, diff_siv
     405    REAL, DIMENSION(klon,nslay) :: lmt_bils
    381406    ! for surface wind stress
    382407    REAL, DIMENSION(klon) :: u0, v0
    … …  
    386411    ! horizontal diffusion and Ekman local vars
    387412    ! dimension = global domain (klon_glo) instead of // subdomains
    388     REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo, dt_ekman_glo
     413    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo,dt_ekman_glo,dt_gm_glo
    389414    ! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop
    390415    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_tmp, dt_ekman_tmp
    … …  
    444469!
    445470!****************************************************************************************
    446     lmt_bils(:)=0.
    447471    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv)
    448472    ! lmt_bils and diff_sst,siv saved by limit_slab
    449     qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400.
    450473    ! qflux = total QFlux correction (in W/m2)
     474    dt_qflux(:,1)=lmt_bils(:,1)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400.
     475    IF (nslay.GT.1) THEN
     476       dt_qflux(:,2:nslay)=lmt_bils(:,2:nslay)
     477    END IF
    451478
    452479!****************************************************************************************
    … …  
    469496        CALL gather(tslab,tslab_glo)
    470497      ! Compute horiz transport on one process only
    471 !$OMP MASTER  ! Only master thread
    472         IF (is_mpi_root) THEN ! Only master processus         
     498        IF (is_mpi_root .AND. is_omp_root) THEN ! Only master processus         
    473499          IF (slab_hdiff) THEN
    474500              dt_hdiff_glo(:,:)=0.
    … …  
    476502          IF (slab_ekman.GT.0) THEN
    477503              dt_ekman_glo(:,:)=0.
     504          END IF
     505          IF (slab_gm) THEN
     506              dt_gm_glo(:,:)=0.
    478507          END IF
    479508          DO i=1,cpl_pas ! time splitting for numerical stability
    … …  
    483512                  CALL slab_ekman1(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
    484513                CASE (2)
    485                   CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
     514                  CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp,dt_hdiff_tmp,slab_gm)
    486515                CASE DEFAULT
    487516                  dt_ekman_tmp(:,:)=0.
    … …  
    493522              ENDDO
    494523              tslab_glo=tslab_glo+dt_ekman_tmp*dtime
     524              IF (slab_gm) THEN ! Gent-McWilliams eddy advection
     525                dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:)
     526                ! convert dt from K.s-1.(kg.m-2) to K.s-1   
     527                DO k=1,nslay
     528                  dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/(slabh(k)*sea_den)
     529                END DO
     530                tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
     531              END IF
    495532            ENDIF
     533! GM included in Ekman_2
     534!            IF (slab_gm) THEN ! Gent-McWilliams eddy advection
     535!              CALL slab_gmdiff(tslab_glo,dt_hdiff_tmp)
     536!              ! convert dt_gm from K.m.s-1 to K.s-1   
     537!              DO k=1,nslay
     538!                dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/slabh(k)
     539!              END DO
     540!              tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
     541!              dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:)
     542!            END IF
    496543            IF (slab_hdiff) THEN ! horizontal diffusion
    497544              ! laplacian of slab T
    … …  
    509556            END DO
    510557          END IF
     558          IF (slab_gm) THEN
     559            !dt_hdiff_glo saved in W/m2
     560            dt_gm_glo(:,:)=dt_gm_glo(:,:)*sea_cap/cpl_pas
     561          END IF
    511562          IF (slab_ekman.GT.0) THEN
    512563            ! dt_ekman_glo saved in W/m2
    513564            dt_ekman_glo(:,:)=dt_ekman_glo(:,:)*sea_cap/cpl_pas
    514565          END IF
    515         END IF ! is_mpi_root
    516 !$OMP END MASTER
     566        END IF ! master process
    517567!$OMP BARRIER
    518568        ! Send new fields back to all processes
    … …  
    520570        IF (slab_hdiff) THEN
    521571          CALL Scatter(dt_hdiff_glo,dt_hdiff)
     572        END IF
     573        IF (slab_gm) THEN
     574          CALL Scatter(dt_gm_glo,dt_gm)
    522575        END IF
    523576        IF (slab_ekman.GT.0) THEN
    … …  
    533586! ***********************************
    534587      ! Add read QFlux
    535       tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas
     588      DO k=1,nslay
     589        tslab(:,k)=tslab(:,k)+dt_qflux(:,k)*cyang*dtime*cpl_pas &
     590                   *slabh(1)/slabh(k)
     591      END DO
    536592      ! Add cumulated surface fluxes
    537593      tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
    … …  
    913969    IF (ALLOCATED(tice)) DEALLOCATE(tice)
    914970    IF (ALLOCATED(seaice)) DEALLOCATE(seaice)
    915     IF (ALLOCATED(slab_bils)) DEALLOCATE(slab_bils)
    916971    IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
    917972    IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
    … …  
    921976    IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman)
    922977    IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff)
     978    IF (ALLOCATED(dt_gm)) DEALLOCATE(dt_gm)
     979    IF (ALLOCATED(dt_qflux)) DEALLOCATE(dt_qflux)
    923980
    924981  END SUBROUTINE ocean_slab_final

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

    r2993 r3002  
    706706  TYPE(ctrl_out), SAVE :: o_tslab = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11/), &
    707707    'tslab', 'Temperature ocean slab', 'K', (/ ('', i=1, 10) /))
     708  TYPE(ctrl_out), SAVE :: o_tslab1 = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11 /), &
     709    'tslab1', 'Temperature ocean slab', 'K', (/ ('', i=1, 10) /))
     710  TYPE(ctrl_out), SAVE :: o_tslab2 = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11 /), &
     711    'tslab2', 'Temperature ocean slab', 'K', (/ ('', i=1, 10) /))
    708712  TYPE(ctrl_out), SAVE :: o_slab_tice = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11/), &
    709713    'slab_tice', 'Temperature banquise slab', 'K', (/ ('', i=1, 10) /))
    … …  
    712716  TYPE(ctrl_out), SAVE :: o_slab_hdiff = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11/), &
    713717    'slab_hdiff', 'Horizontal diffusion', 'W/m2', (/ ('', i=1, 10) /))
     718  TYPE(ctrl_out), SAVE :: o_slab_gm = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11 /), &
     719    'slab_gm', 'GM eddy advection', 'W/m2', (/ ('', i=1, 10) /))
    714720  TYPE(ctrl_out), SAVE :: o_slab_ekman = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11/), &
    715721    'slab_ekman', 'Ekman heat transport', 'W/m2', (/ ('', i=1, 10) /))

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

    r2993 r3002  
    9191         o_slab_qflux, o_tslab, o_slab_bils, &
    9292         o_slab_bilg, o_slab_sic, o_slab_tice, &
    93          o_slab_hdiff, o_slab_ekman, &
     93         o_slab_hdiff, o_slab_ekman, o_slab_gm,  &
    9494         o_weakinv, o_dthmin, o_cldtau, &
    9595         o_cldemi, o_pr_con_l, o_pr_con_i, &
    … …  
    329329         ok_4xCO2atm
    330330
    331     USE ocean_slab_mod, ONLY: nslay, tslab, slab_bils, slab_bilg, tice, &
    332         seaice, slab_ekman,slab_hdiff, dt_ekman, dt_hdiff
     331    USE ocean_slab_mod, ONLY: nslay, tslab, slab_bilg, tice, seaice, &
     332        slab_ekman,slab_hdiff,slab_gm,dt_ekman, dt_hdiff, dt_gm, dt_qflux
    333333    USE pbl_surface_mod, ONLY: snow
    334334    USE indice_sol_mod, ONLY: nbsrf
    … …  
    11031103       ! Output of slab ocean variables
    11041104       IF (type_ocean=='slab ') THEN
    1105           CALL histwrite_phy(o_slab_qflux, slab_wfbils)
    1106           CALL histwrite_phy(o_slab_bils, slab_bils)
     1105          CALL histwrite_phy(o_slab_bils, slab_wfbils)
    11071106          IF (nslay.EQ.1) THEN
    11081107              zx_tmp_fi2d(:)=tslab(:,1)
    11091108              CALL histwrite_phy(o_tslab, zx_tmp_fi2d)
     1109              zx_tmp_fi2d(:)=dt_qflux(:,1)
     1110              CALL histwrite_phy(o_slab_qflux, zx_tmp_fi2d)
    11101111          ELSE
    11111112              CALL histwrite_phy(o_tslab, tslab(:,1:nslay))
     1113              CALL histwrite_phy(o_slab_qflux, dt_qflux(:,1:nslay))
    11121114          ENDIF
    11131115          IF (version_ocean=='sicINT') THEN
    … …  
    11161118              CALL histwrite_phy(o_slab_sic, seaice)
    11171119          ENDIF
     1120          IF (slab_gm) THEN
     1121             CALL histwrite_phy(o_slab_gm, dt_gm(:,1:nslay))
     1122          END IF
    11181123          IF (slab_hdiff) THEN
    11191124            IF (nslay.EQ.1) THEN

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

    r2656 r3002  
    1111  REAL,SAVE,ALLOCATABLE :: fext(:) ! Coriolis f times cell area
    1212  !$OMP THREADPRIVATE(fext)
     13  REAL,SAVE,ALLOCATABLE :: beta(:) ! df/dy
     14  !$OMP THREADPRIVATE(beta)
    1315  REAL,SAVE,ALLOCATABLE :: unsairez(:) ! 1/cell area
    1416  !$OMP THREADPRIVATE(unsairez)
    … …  
    3436  !$OMP THREADPRIVATE(airev)
    3537
    36   ! Local variables for horiz mass flux in slab
    37   LOGICAL,SAVE :: alpha_var
     38  ! Local parameters for slab transport
     39  LOGICAL,SAVE :: alpha_var ! variable coef for deep temp (1 layer)
    3840  !$OMP THREADPRIVATE(alpha_var)
    39   LOGICAL,SAVE :: slab_upstream
     41  LOGICAL,SAVE :: slab_upstream ! upstream scheme ? (1 layer)
    4042  !$OMP THREADPRIVATE(slab_upstream)
    41 
     43  LOGICAL,SAVE :: slab_sverdrup ! use wind stress curl at equator
     44  !$OMP THREADPRIVATE(slab_sverdrup)
     45  LOGICAL,SAVE :: slab_tyeq ! use merid wind stress at equator
     46  !$OMP THREADPRIVATE(slab_tyeq)
     47  LOGICAL,SAVE :: ekman_zonadv ! use zonal advection by Ekman currents
     48  !$OMP THREADPRIVATE(ekman_zonadv)
     49  LOGICAL,SAVE :: ekman_zonavg ! zonally average wind stress
     50  !$OMP THREADPRIVATE(ekman_zonavg)
     51
     52  REAL,SAVE :: alpham
     53  !$OMP THREADPRIVATE(alpham)
     54  REAL,SAVE :: gmkappa
     55  !$OMP THREADPRIVATE(gmkappa)
     56  REAL,SAVE :: gm_smax
     57  !$OMP THREADPRIVATE(gm_smax)
     58
     59! geometry variables : f, beta, mask...
    4260  REAL,SAVE,ALLOCATABLE :: zmasqu(:) ! continent mask for zonal mass flux
    4361  !$OMP THREADPRIVATE(zmasqu)
    … …  
    4664  REAL,SAVE,ALLOCATABLE :: unsfv(:) ! 1/f, v points
    4765  !$OMP THREADPRIVATE(unsfv)
     66  REAL,SAVE,ALLOCATABLE :: unsbv(:) ! 1/beta
     67  !$OMP THREADPRIVATE(unsbv)
    4868  REAL,SAVE,ALLOCATABLE :: unsev(:) ! 1/epsilon (drag)
    4969  !$OMP THREADPRIVATE(unsev)
    … …  
    6383                                  cu_,cuvsurcv_,cv_,cvusurcu_, &
    6484                                  aire_,apoln_,apols_, &
    65                                   aireu_,airev_)
     85                                  aireu_,airev_,rlatv)
     86    USE comconst_mod, ONLY: omeg, rad
    6687    ! number of points in lon, lat
    6788    IMPLICIT NONE
    … …  
    82103    REAL,INTENT(IN) :: aireu_(ip1jmp1)
    83104    REAL,INTENT(IN) :: airev_(ip1jm)
    84 
     105    REAL,INTENT(IN) :: rlatv(nbp_lat-1)
    85106
    86107    ! Sanity check on dimensions
    … …  
    112133    aireu(:)=aireu_(:)
    113134    allocate(airev(ip1jm))
    114     airev(:)=airev_(:)
     135    airev(:)=airev_(:)
     136    allocate(beta(nbp_lat-1))
     137    beta(:)=2*omeg*cos(rlatv(:))/rad
    115138
    116139  END SUBROUTINE ini_slab_transp_geom
    … …  
    139162    ip1jmp1=(nbp_lon+1)*nbp_lat
    140163
     164! Options for Heat transport
     165    ! Alpha variable?
     166      alpha_var=.FALSE.
     167      CALL getin('slab_alphav',alpha_var)
     168      print *,'alpha variable',alpha_var
     169!  centered ou upstream scheme for meridional transport
     170      slab_upstream=.FALSE.
     171      CALL getin('slab_upstream',slab_upstream)
     172      print *,'upstream slab scheme',slab_upstream
     173! Sverdrup balance at equator ?
     174      slab_sverdrup=.TRUE.
     175      CALL getin('slab_sverdrup',slab_sverdrup)
     176      print *,'Sverdrup balance',slab_sverdrup
     177! Use tauy for meridional flux at equator ?
     178      slab_tyeq=.TRUE.
     179      CALL getin('slab_tyeq',slab_tyeq)
     180      print *,'Tauy forcing at equator',slab_tyeq
     181! Use tauy for meridional flux at equator ?
     182      ekman_zonadv=.TRUE.
     183      CALL getin('slab_ekman_zonadv',ekman_zonadv)
     184      print *,'Use Ekman flow in zonal direction',ekman_zonadv
     185! Use tauy for meridional flux at equator ?
     186      ekman_zonavg=.FALSE.
     187      CALL getin('slab_ekman_zonavg',ekman_zonavg)
     188      print *,'Use zonally-averaged wind stress ?',ekman_zonavg
     189! Value of alpha
     190      alpham=2./3.
     191      CALL getin('slab_alpha',alpham)
     192      print *,'slab_alpha',alpham
     193! GM k coefficient (m2/s) for 2-layers
     194      gmkappa=1000.
     195      CALL getin('slab_gmkappa',gmkappa)
     196      print *,'slab_gmkappa',gmkappa
     197! GM k coefficient (m2/s) for 2-layers
     198      gm_smax=2e-3
     199      CALL getin('slab_gm_smax',gm_smax)
     200      print *,'slab_gm_smax',gm_smax
     201! -----------------------------------------------------------
    141202! Define ocean / continent mask (no flux into continent cell)
    142203    allocate(zmasqu(ip1jmp1))
    … …  
    161222      END IF
    162223    END DO
     224
     225! -----------------------------------------------------------
     226! Coriolis and friction for Ekman transport
    163227    slab_ekman=2
    164228    CALL getin("slab_ekman",slab_ekman)
    165 ! Coriolis and friction for Ekman transport
    166229    IF (slab_ekman.GT.0) THEN
    167230      allocate(unsfv(ip1jm))
    … …  
    169232      allocate(unsfu(ip1jmp1))
    170233      allocate(unseu(ip1jmp1))
     234      allocate(unsbv(ip1jm))
    171235
    172236      eps=1e-5 ! Drag
    … …  
    176240      ! coefs to convert tau_x, tau_y to Ekman mass fluxes
    177241      ! on 2D grid v points
     242      ! Compute correction factor [0 1] near the equator (f<<eps)
     243      IF (slab_sverdrup) THEN
     244         ! New formulation, sharper near equator, when eps gives Rossby Radius
     245         DO i=1,ip1jm
     246           unsev(i)=exp(-ff(i)*ff(i)/eps**2)
     247         ENDDO
     248      ELSE
     249         DO i=1,ip1jm
     250           unsev(i)=eps**2/(ff(i)*ff(i)+eps**2)
     251         ENDDO
     252      END IF ! slab_sverdrup
     253      ! 1/beta
     254      DO i=1,jjp1-1
     255        unsbv((i-1)*iip1+1:i*iip1)=unsev((i-1)*iip1+1:i*iip1)/beta(i)
     256      END DO
     257      ! 1/f
     258      ff=SIGN(MAX(ABS(ff),eps/100.),ff) ! avoid value 0 at equator...
    178259      DO i=1,ip1jm
    179         unsev(i)=eps/(ff(i)*ff(i)+eps**2)
    180         unsfv(i)=ff(i)/(ff(i)*ff(i)+eps**2)
    181       ENDDO
     260        unsfv(i)=(1.-unsev(i))/ff(i)
     261      END DO
    182262      ! compute values on 2D u grid
     263      ! 1/eps
     264      unsev(:)=unsev(:)/eps
    183265      CALL gr_v_scal(1,unsfv,unsfu)
    184266      CALL gr_v_scal(1,unsev,unseu)
    185267    END IF
    186 
    187 ! Other options for Ekman transport
    188     ! Alpha variable?
    189       alpha_var=.FALSE.
    190       CALL getin('slab_alphav',alpha_var)
    191       print *,'alpha variable',alpha_var
    192 !  centered ou upstream scheme for meridional transport
    193       slab_upstream=.FALSE.
    194       CALL getin('slab_upstream',slab_upstream)
    195       print *,'upstream slab scheme',slab_upstream
    196268 
    197269  END SUBROUTINE ini_slab_transp
    … …  
    249321      REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat)
    250322      REAL txv((nbp_lon+1)*(nbp_lat-1)),tyv((nbp_lon+1)*(nbp_lat-1))
     323      REAL tcurl((nbp_lon+1)*(nbp_lat-1))
    251324      ! zonal and meridional Ekman mass fluxes at u, v points (2D grid)
    252325      REAL fluxz((nbp_lon+1)*nbp_lat),fluxm((nbp_lon+1)*(nbp_lat-1))
    … …  
    272345! First and last points in zonal direction are the same
    273346! we use 1 index ij from 1 to (iim+1)*(jjm+1)
     347      ! north and south poles
     348      tx_phy(1)=0.
     349      tx_phy(klon_glo)=0.
     350      ty_phy(1)=0.
     351      ty_phy(klon_glo)=0.
    274352      CALL gr_fi_dyn(1,iip1,jjp1,tx_phy,txu)
    275353      CALL gr_fi_dyn(1,iip1,jjp1,ty_phy,tyu)
    … …  
    278356      ! Meridional mass flux
    279357      CALL gr_scal_v(1,txu,txv) ! wind stress at v points
    280       CALL gr_scal_v(1,tyu,tyv)
    281       fluxm=tyv*unsev-txv*unsfv ! in kg.s-1.m-1 (zonal distance)
     358      IF (slab_sverdrup) THEN ! Sverdrup bal. near equator
     359        tcurl=(txu(1:ip1jm)-txu(iip2:ip1jmp1))/cv(:)
     360        fluxm=-tcurl*unsbv-txv*unsfv ! in kg.s-1.m-1 (zonal distance)
     361      ELSE
     362        CALL gr_scal_v(1,tyu,tyv)
     363        fluxm=tyv*unsev-txv*unsfv ! in kg.s-1.m-1 (zonal distance)
     364      ENDIF
    282365      ! Zonal mass flux
    283366      CALL gr_scal_u(1,txu,txu) ! wind stress at u points
    … …  
    316399          ! to avoid "hot spots" where there is surface convergence
    317400          DO ij=iip2,ip1jm
    318               alpha(ij)=2./3.-fluxv(ij)/fluxt(ij)/3.
     401              alpha(ij)=alpham-fluxv(ij)/fluxt(ij)*(1.-alpham)
    319402          ENDDO
    320           alpha(1:iip1)=2./3.-fluxv(1)/fluxt(1)/3.
    321           alpha(ip1jm+1:ip1jmp1)=2./3.-fluxv(ip1jmp1)/fluxt(ip1jmp1)/3.
     403          alpha(1:iip1)=alpham-fluxv(1)/fluxt(1)*(1.-alpham)
     404          alpha(ip1jm+1:ip1jmp1)=alpham-fluxv(ip1jmp1)/fluxt(ip1jmp1)*(1.-alpham)
    322405      ELSE
    323           alpha(:)=2./3.
     406          alpha(:)=alpham
    324407          ! Tsurf-Tdeep ~ 10° in the Tropics
    325408      ENDIF
    … …  
    380463  END SUBROUTINE slab_ekman1
    381464
    382   SUBROUTINE slab_ekman2(tx_phy,ty_phy,ts_phy,dt_phy)
     465  SUBROUTINE slab_ekman2(tx_phy,ty_phy,ts_phy,dt_phy_ek,dt_phy_gm,slab_gm)
    383466! Temperature tendency for 2-layers slab ocean
    384467! note : tendency dt later multiplied by (delta time)/(rho.H)
    385468! to convert from divergence of heat fluxes to T
     469
     470! 11/16 : Inclusion of GM-like eddy advection
     471
     472      IMPLICIT NONE
     473     
     474      LOGICAL,INTENT(in) :: slab_gm
     475      ! Here, temperature and flux variables are on 2 layers
     476      INTEGER ij
     477
     478      ! wind stress variables
     479      REAL tx_phy(klon_glo),ty_phy(klon_glo)
     480      REAL txv((nbp_lon+1)*(nbp_lat-1)), tyv((nbp_lon+1)*(nbp_lat-1))
     481      REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat)
     482      REAL tcurl((nbp_lon+1)*(nbp_lat-1))
     483      ! slab temperature on  1D, 2D grid
     484      REAL ts_phy(klon_glo,2), ts((nbp_lon+1)*nbp_lat,2)
     485      ! Temperature gradient, v-points
     486      REAL dty((nbp_lon+1)*(nbp_lat-1)),dtx((nbp_lon+1)*nbp_lat)
     487      ! Vertical temperature difference, V-points
     488      REAL dtz((nbp_lon+1)*(nbp_lat-1))
     489      ! zonal and meridional mass fluxes at u, v points (2D grid)
     490      REAL fluxz((nbp_lon+1)*nbp_lat), fluxm((nbp_lon+1)*(nbp_lat-1))
     491      ! vertical mass flux between the 2 layers
     492      REAL fluxv_ek((nbp_lon+1)*nbp_lat)
     493      REAL fluxv_gm((nbp_lon+1)*nbp_lat)
     494      ! zonal and meridional heat fluxes
     495      REAL fluxtz((nbp_lon+1)*nbp_lat,2)
     496      REAL fluxtm((nbp_lon+1)*(nbp_lat-1),2)
     497      ! temperature tendency (in K.s-1.kg.m-2)
     498      REAL dt_ek((nbp_lon+1)*nbp_lat,2), dt_phy_ek(klon_glo,2)
     499      REAL dt_gm((nbp_lon+1)*nbp_lat,2), dt_phy_gm(klon_glo,2)
     500      ! helper vars
     501      REAL zonavg, fluxv
     502      REAL, PARAMETER :: sea_den=1025. ! sea water density
     503
     504      INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1
     505
     506! Grid definitions
     507      iim=nbp_lon
     508      iip1=nbp_lon+1
     509      iip2=nbp_lon+2
     510      jjp1=nbp_lat
     511      ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
     512      ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
     513      ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
     514! Convert temperature to 2D grid
     515      CALL gr_fi_dyn(2,iip1,jjp1,ts_phy,ts)
     516
     517! ------------------------------------
     518! Ekman mass fluxes and Temp tendency
     519! ------------------------------------
     520! Convert taux,y to 2D  scalar grid
     521      ! north and south poles tx,ty no meaning
     522      tx_phy(1)=0.
     523      tx_phy(klon_glo)=0.
     524      ty_phy(1)=0.
     525      ty_phy(klon_glo)=0.
     526      CALL gr_fi_dyn(1,iip1,jjp1,tx_phy,txu)
     527      CALL gr_fi_dyn(1,iip1,jjp1,ty_phy,tyu)
     528      IF (ekman_zonavg) THEN ! use zonal average of wind stress
     529        DO ij=1,jjp1-2
     530          zonavg=SUM(txu(ij*iip1+1:ij*iip1+iim))/iim
     531          txu(ij*iip1+1:(ij+1)*iip1)=zonavg
     532          zonavg=SUM(tyu(ij*iip1+1:ij*iip1+iim))/iim
     533          tyu(ij*iip1+1:(ij+1)*iip1)=zonavg
     534        END DO
     535      END IF
     536         
     537! Divide taux,y by f or eps, and convert to 2D u,v grids
     538! (Arakawa C grid)
     539      ! Meridional flux
     540      CALL gr_scal_v(1,txu,txv) ! wind stress at v points
     541      fluxm=-txv*unsfv ! in kg.s-1.m-1 (zonal distance)
     542      IF (slab_sverdrup) THEN ! Sverdrup bal. near equator
     543        tcurl=(txu(1:ip1jm)-txu(iip2:ip1jmp1))/cv(:) ! dtx/dy
     544        !poles curl = 0
     545        tcurl(1:iip1)=0.
     546        tcurl(ip1jmi1+1:ip1jm)=0.
     547        fluxm=fluxm-tcurl*unsbv
     548      ENDIF
     549      IF (slab_tyeq) THEN ! meridional wind forcing at equator
     550        CALL gr_scal_v(1,tyu,tyv)
     551        fluxm=fluxm+tyv*unsev ! in kg.s-1.m-1 (zonal distance)
     552      ENDIF
     553!  apply continent mask, multiply by horiz grid dimension
     554      fluxm=fluxm*cv*cuvsurcv*zmasqv
     555
     556      ! Zonal flux
     557      IF (ekman_zonadv) THEN
     558        CALL gr_scal_u(1,txu,txu) ! wind stress at u points
     559        CALL gr_scal_u(1,tyu,tyu)
     560        fluxz=tyu*unsfu+txu*unseu
     561        !  apply continent mask, multiply by horiz grid dimension
     562        fluxz=fluxz*cu*cvusurcu*zmasqu
     563      END IF
     564           
     565!  Vertical mass flux from mass budget (divergence of horiz fluxes)
     566      IF (ekman_zonadv) THEN
     567         DO ij=iip2,ip1jm
     568           fluxv_ek(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1)
     569         ENDDO
     570      ELSE
     571         DO ij=iip2,ip1jm
     572           fluxv_ek(ij)=-fluxm(ij)+fluxm(ij-iip1)
     573         ENDDO
     574      END IF
     575      DO ij=iip1,ip1jmi1,iip1
     576         fluxv_ek(ij+1)=fluxv_ek(ij+iip1)
     577      END DO
     578!  vertical mass flux at Poles
     579      fluxv_ek(1)=-SUM(fluxm(1:iim))     
     580      fluxv_ek(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1))
     581
     582! Meridional heat fluxes
     583      DO ij=1,ip1jm
     584          ! centered scheme
     585          fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2.
     586          fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2.
     587      END DO
     588
     589! Zonal heat fluxes
     590! Schema upstream     
     591      IF (ekman_zonadv) THEN
     592        DO ij=iip2,ip1jm
     593          IF (fluxz(ij).GE.0.) THEN
     594                 fluxtz(ij,1)=fluxz(ij)*ts(ij,1)
     595                 fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2)
     596          ELSE
     597                 fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1)
     598                 fluxtz(ij,2)=-fluxz(ij)*ts(ij,2)
     599          ENDIF
     600        ENDDO
     601        DO ij=iip1*2,ip1jmp1,iip1
     602               fluxtz(ij,:)=fluxtz(ij-iim,:)
     603        END DO
     604      ELSE
     605        fluxtz(:,:)=0.
     606      ENDIF       
     607                   
     608! Temperature tendency, horizontal advection:
     609      DO ij=iip2,ip1jm
     610         dt_ek(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) &
     611                 +fluxtm(ij,:)-fluxtm(ij-iip1,:)
     612      END DO
     613! Poles
     614      dt_ek(1,:)=SUM(fluxtm(1:iim,:),dim=1)
     615      dt_ek(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1)
     616
     617! ------------------------------------
     618! GM mass fluxes and Temp tendency
     619! ------------------------------------
     620     IF (slab_gm) THEN
     621! Vertical Temperature difference T1-T2 on v-grid points
     622      CALL gr_scal_v(1,ts(:,1)-ts(:,2),dtz)
     623      dtz(:)=MAX(dtz(:),0.25)
     624! Horizontal Temperature differences
     625      CALL grad(1,(ts(:,1)+ts(:,2))/2.,dtx,dty)
     626! Meridional flux = -k.s (s=dyT/dzT)
     627! Continent mask, multiply by dz/dy
     628      fluxm=dty/dtz*500.*cuvsurcv*zmasqv
     629! slope limitation, multiply by kappa
     630      fluxm=-gmkappa*SIGN(MIN(ABS(fluxm),gm_smax*cv*cuvsurcv),dty)
     631! convert to kg/s
     632      fluxm(:)=fluxm(:)*sea_den
     633! Zonal flux = 0. (temporary)
     634      fluxz(:)=0.
     635!  Vertical mass flux from mass budget (divergence of horiz fluxes)
     636      DO ij=iip2,ip1jm
     637        fluxv_gm(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1)
     638      ENDDO
     639      DO ij=iip1,ip1jmi1,iip1
     640         fluxv_gm(ij+1)=fluxv_gm(ij+iip1)
     641      END DO
     642!  vertical mass flux at Poles
     643      fluxv_gm(1)=-SUM(fluxm(1:iim))     
     644      fluxv_gm(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1))
     645
     646! Meridional heat fluxes
     647      DO ij=1,ip1jm
     648          ! centered scheme
     649          fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2.
     650          fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2.
     651      END DO
     652
     653! Zonal heat fluxes
     654! Schema upstream     
     655      DO ij=iip2,ip1jm
     656      IF (fluxz(ij).GE.0.) THEN
     657             fluxtz(ij,1)=fluxz(ij)*ts(ij,1)
     658             fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2)
     659      ELSE
     660             fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1)
     661             fluxtz(ij,2)=-fluxz(ij)*ts(ij,2)
     662      ENDIF
     663      ENDDO
     664      DO ij=iip1*2,ip1jmp1,iip1
     665             fluxtz(ij,:)=fluxtz(ij-iim,:)
     666      END DO
     667                   
     668! Temperature tendency :
     669! divergence of horizontal heat fluxes
     670      DO ij=iip2,ip1jm
     671         dt_gm(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) &
     672                 +fluxtm(ij,:)-fluxtm(ij-iip1,:)
     673      END DO
     674! Poles
     675      dt_gm(1,:)=SUM(fluxtm(1:iim,:),dim=1)
     676      dt_gm(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1)
     677     ELSE
     678       dt_gm(:,:)=0.
     679       fluxv_gm(:)=0.
     680     ENDIF ! slab_gm
     681
     682! ------------------------------------
     683! Temp tendency from vertical advection
     684! Divide by cell area
     685! ------------------------------------
     686! vertical heat flux = mass flux * T, upstream scheme
     687      DO ij=iip2,ip1jm
     688         fluxv=fluxv_ek(ij)+fluxv_gm(ij) ! net flux, needed for upstream scheme
     689         IF (fluxv.GT.0.) THEN
     690           dt_ek(ij,1)=dt_ek(ij,1)+fluxv_ek(ij)*ts(ij,2)
     691           dt_ek(ij,2)=dt_ek(ij,2)-fluxv_ek(ij)*ts(ij,2)
     692           dt_gm(ij,1)=dt_gm(ij,1)+fluxv_gm(ij)*ts(ij,2)
     693           dt_gm(ij,2)=dt_gm(ij,2)-fluxv_gm(ij)*ts(ij,2)
     694         ELSE
     695           dt_ek(ij,1)=dt_ek(ij,1)+fluxv_ek(ij)*ts(ij,1)
     696           dt_ek(ij,2)=dt_ek(ij,2)-fluxv_ek(ij)*ts(ij,1)
     697           dt_gm(ij,1)=dt_gm(ij,1)+fluxv_gm(ij)*ts(ij,1)
     698           dt_gm(ij,2)=dt_gm(ij,2)-fluxv_gm(ij)*ts(ij,1)
     699         ENDIF
     700         ! divide by cell area
     701         dt_ek(ij,:)=dt_ek(ij,:)/aire(ij)
     702         dt_gm(ij,:)=dt_gm(ij,:)/aire(ij)
     703      END DO
     704      ! North Pole
     705      fluxv=fluxv_ek(1)+fluxv_gm(1)
     706        IF (fluxv.GT.0.) THEN
     707          dt_ek(1,1)=dt_ek(1,1)+fluxv_ek(1)*ts(1,2)
     708          dt_ek(1,2)=dt_ek(1,2)-fluxv_ek(1)*ts(1,2)
     709          dt_gm(1,1)=dt_gm(1,1)+fluxv_gm(1)*ts(1,2)
     710          dt_gm(1,2)=dt_gm(1,2)-fluxv_gm(1)*ts(1,2)
     711        ELSE
     712          dt_ek(1,1)=dt_ek(1,1)+fluxv_ek(1)*ts(1,1)
     713          dt_ek(1,2)=dt_ek(1,2)-fluxv_ek(1)*ts(1,1)
     714          dt_gm(1,1)=dt_gm(1,1)+fluxv_gm(1)*ts(1,1)
     715          dt_gm(1,2)=dt_gm(1,2)-fluxv_gm(1)*ts(1,1)
     716        ENDIF
     717      dt_ek(1,:)=dt_ek(1,:)/apoln
     718      dt_gm(1,:)=dt_gm(1,:)/apoln
     719      ! South pole
     720      fluxv=fluxv_ek(ip1jmp1)+fluxv_gm(ip1jmp1)
     721        IF (fluxv.GT.0.) THEN
     722          dt_ek(ip1jmp1,1)=dt_ek(ip1jmp1,1)+fluxv_ek(ip1jmp1)*ts(ip1jmp1,2)
     723          dt_ek(ip1jmp1,2)=dt_ek(ip1jmp1,2)-fluxv_ek(ip1jmp1)*ts(ip1jmp1,2)
     724          dt_gm(ip1jmp1,1)=dt_gm(ip1jmp1,1)+fluxv_gm(ip1jmp1)*ts(ip1jmp1,2)
     725          dt_gm(ip1jmp1,2)=dt_gm(ip1jmp1,2)-fluxv_gm(ip1jmp1)*ts(ip1jmp1,2)
     726        ELSE
     727          dt_ek(ip1jmp1,1)=dt_ek(ip1jmp1,1)+fluxv_ek(ip1jmp1)*ts(ip1jmp1,1)
     728          dt_ek(ip1jmp1,2)=dt_ek(ip1jmp1,2)-fluxv_ek(ip1jmp1)*ts(ip1jmp1,1)
     729          dt_gm(ip1jmp1,1)=dt_gm(ip1jmp1,1)+fluxv_gm(ip1jmp1)*ts(ip1jmp1,1)
     730          dt_gm(ip1jmp1,2)=dt_gm(ip1jmp1,2)-fluxv_gm(ip1jmp1)*ts(ip1jmp1,1)
     731        ENDIF
     732      dt_ek(ip1jmp1,:)=dt_ek(ip1jmp1,:)/apols
     733      dt_gm(ip1jmp1,:)=dt_gm(ip1jmp1,:)/apols
     734     
     735      dt_ek(2:iip1,1)=dt_ek(1,1)
     736      dt_ek(2:iip1,2)=dt_ek(1,2)
     737      dt_gm(2:iip1,1)=dt_gm(1,1)
     738      dt_gm(2:iip1,2)=dt_gm(1,2)
     739      dt_ek(ip1jm+1:ip1jmp1,1)=dt_ek(ip1jmp1,1)
     740      dt_ek(ip1jm+1:ip1jmp1,2)=dt_ek(ip1jmp1,2)
     741      dt_gm(ip1jm+1:ip1jmp1,1)=dt_gm(ip1jmp1,1)
     742      dt_gm(ip1jm+1:ip1jmp1,2)=dt_gm(ip1jmp1,2)
     743
     744      DO ij=iip1,ip1jmi1,iip1
     745         dt_gm(ij+1,:)=dt_gm(ij+iip1,:)
     746         dt_ek(ij+1,:)=dt_ek(ij+iip1,:)
     747      END DO
     748
     749! T tendency back to 1D grid...
     750      CALL gr_dyn_fi(2,iip1,jjp1,dt_ek,dt_phy_ek)
     751      CALL gr_dyn_fi(2,iip1,jjp1,dt_gm,dt_phy_gm)
     752
     753      RETURN
     754  END SUBROUTINE slab_ekman2
     755
     756  SUBROUTINE slab_gmdiff(ts_phy,dt_phy)
     757! Temperature tendency for 2-layers slab ocean
     758! Due to Gent-McWilliams type eddy-induced advection
    386759     
    387760      IMPLICIT NONE
    … …  
    389762      ! Here, temperature and flux variables are on 2 layers
    390763      INTEGER ij
    391 
    392       REAL tx_phy(klon_glo),ty_phy(klon_glo)
    393       REAL txv((nbp_lon+1)*(nbp_lat-1)), tyv((nbp_lon+1)*(nbp_lat-1))
    394       REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat)
     764      ! Temperature gradient, v-points
     765      REAL dty((nbp_lon+1)*(nbp_lat-1)),dtx((nbp_lon+1)*nbp_lat)
     766      ! Vertical temperature difference, V-points
     767      REAL dtz((nbp_lon+1)*(nbp_lat-1))
    395768      ! slab temperature on  1D, 2D grid
    396       REAL ts_phy(klon_glo,2), ts((nbp_lon+1)*nbp_lat,2)
    397       ! zonal and meridional Ekman mass flux at u, v points (2D grid)
     769      REAL ts_phy(klon_glo,2),ts((nbp_lon+1)*nbp_lat,2)
     770      ! zonal and meridional mass fluxes at u, v points (2D grid)
    398771      REAL fluxz((nbp_lon+1)*nbp_lat), fluxm((nbp_lon+1)*(nbp_lat-1))
    399772      ! vertical mass flux between the 2 layers
    … …  
    416789      ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
    417790
    418 ! Convert taux,y to 2D  scalar grid
    419       CALL gr_fi_dyn(1,iip1,jjp1,tx_phy,txu)
    420       CALL gr_fi_dyn(1,iip1,jjp1,ty_phy,tyu)
    421 ! Multiply taux,y by f or eps, and convert to 2D u,v grids
    422 ! (Arakawa C grid)
    423       ! Meridional flux
    424       CALL gr_scal_v(1,txu,txv) ! wind stress at v points
    425       CALL gr_scal_v(1,tyu,tyv)
    426       fluxm=tyv*unsev-txv*unsfv
    427       ! Zonal flux
    428       CALL gr_scal_u(1,txu,txu) ! wind stress at u points
    429       CALL gr_scal_u(1,tyu,tyu)
    430       fluxz=tyu*unsfu+txu*unseu
    431            
    432791! Convert temperature to 2D grid
    433792      CALL gr_fi_dyn(2,iip1,jjp1,ts_phy,ts)
    434        
    435 !  Mass fluxes (apply continent mask, multiply by horiz grid dimension)
    436       fluxm=fluxm*cv*cuvsurcv*zmasqv
    437       fluxz=fluxz*cu*cvusurcu*zmasqu
     793! Vertical Temperature difference T1-T2 on v-grid points
     794      CALL gr_scal_v(1,ts(:,1)-ts(:,2),dtz)
     795      dtz(:)=MAX(dtz(:),0.25)
     796! Horizontal Temperature differences
     797      CALL grad(1,(ts(:,1)+ts(:,2))/2.,dtx,dty)
     798! Meridional flux = -k.s (s=dyT/dzT)
     799! Continent mask, multiply by dz/dy
     800      fluxm=dty/dtz*500.*cuvsurcv*zmasqv
     801! slope limitation, multiply by kappa
     802      fluxm=-gmkappa*SIGN(MIN(ABS(fluxm),gm_smax*cv*cuvsurcv),dty)
     803! Zonal flux = 0. (temporary)
     804      fluxz(:)=0.
    438805!  Vertical mass flux from mass budget (divergence of horiz fluxes)
    439806      DO ij=iip2,ip1jm
    … …  
    489856         dt(ij+1,:)=dt(ij+iip1,:)
    490857      END DO
    491 ! Pôles
     858! Poles
    492859      dt(1,:)=SUM(fluxtm(1:iim,:),dim=1)
    493860        IF (fluxv(1).GT.0.) THEN
    … …  
    517884
    518885      RETURN
    519   END SUBROUTINE slab_ekman2
     886  END SUBROUTINE slab_gmdiff
    520887
    521888!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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

    r2719 r3002  
    2929  USE ocean_cpl_mod, ONLY    : ocean_cpl_noice
    3030  USE indice_sol_mod, ONLY : nbsrf, is_oce
     31  USE limit_read_mod
    3132!
    3233! This subroutine will make a call to ocean_XXX_noice according to the ocean mode (force,
    … …  
    9192    REAL                  :: tmp
    9293    REAL, PARAMETER       :: cepdu2=(0.1)**2
    93     REAL, DIMENSION(klon) :: alb_eau
     94    REAL, DIMENSION(klon) :: alb_eau, z0_lim
    9495    REAL, DIMENSION(klon) :: radsol
    9596    REAL, DIMENSION(klon) :: cdragq ! Cdrag pour l'evaporation
    … …  
    123124       cdragq(1:knon)=cdragh(1:knon)
    124125    ENDIF
     126
    125127
    126128!******************************************************************************
    … …  
    224226    alb_dif_new(1:knon,:)=MIN(MAX(alb_dif_new(1:knon,:),0.0),1.0)
    225227!
     228! F. Codron albedo read from limit.nc
     229ELSE IF (iflag_albedo==2) THEN
     230    CALL limit_read_rug_alb(itime, dtime, jour,&
     231         knon, knindex, z0_lim, alb_eau)
     232    DO i =1, knon
     233      DO  k=1,nsw
     234       alb_dir_new(i,k) = alb_eau(i)
     235      ENDDO
     236    ENDDO
     237    alb_dif_new=alb_dir_new
    226238ENDIF
    227239!albedo SB <<<
Note: See TracChangeset for help on using the changeset viewer.