Index: LMDZ5/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90 =================================================================== --- LMDZ5/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90 (revision 3001) +++ LMDZ5/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90 (revision 3002) @@ -131,5 +131,5 @@ cu,cuvsurcv,cv,cvusurcu, & aire,apoln,apols, & - aireu,airev) + aireu,airev,rlatvdyn) END IF Index: LMDZ5/trunk/libf/phylmd/conf_phys_m.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/conf_phys_m.F90 (revision 3001) +++ LMDZ5/trunk/libf/phylmd/conf_phys_m.F90 (revision 3002) @@ -2419,5 +2419,5 @@ !--test on ocean surface albedo - IF (iflag_albedo.LT.0.OR.iflag_albedo.GT.1) THEN + IF (iflag_albedo.LT.0.OR.iflag_albedo.GT.2) THEN WRITE(lunout,*) ' ERROR iflag_albedo<>0,1' CALL abort_physic('conf_phys','choice iflag_albedo not valid',1) Index: LMDZ5/trunk/libf/phylmd/limit_slab.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/limit_slab.F90 (revision 3001) +++ LMDZ5/trunk/libf/phylmd/limit_slab.F90 (revision 3002) @@ -8,4 +8,5 @@ USE netcdf USE indice_sol_mod + USE ocean_slab_mod, ONLY: nslay IMPLICIT NONE @@ -18,10 +19,11 @@ INTEGER, INTENT(IN) :: jour ! jour a lire dans l'annee REAL , INTENT(IN) :: dtime ! pas de temps de la physique (en s) - REAL, DIMENSION(klon), INTENT(OUT) :: lmt_bils, diff_sst, diff_siv + REAL, DIMENSION(klon), INTENT(OUT) :: diff_sst, diff_siv + REAL, DIMENSION(klon,nslay), INTENT(OUT) :: lmt_bils ! Locals variables with attribute SAVE !**************************************************************************************** - REAL, DIMENSION(:), ALLOCATABLE, SAVE :: bils_save, diff_sst_save - REAL, DIMENSION(:), ALLOCATABLE, SAVE :: diff_siv_save + REAL, DIMENSION(:), ALLOCATABLE, SAVE :: diff_siv_save, diff_sst_save + REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: bils_save !$OMP THREADPRIVATE(bils_save, diff_sst_save, diff_siv_save) @@ -31,7 +33,9 @@ INTEGER :: nvarid, nid, ierr, i INTEGER, DIMENSION(2) :: start, epais - REAL, DIMENSION(klon_glo):: bils_glo, sst_l_glo, sst_lp1_glo, diff_sst_glo + REAL, DIMENSION(klon_glo):: sst_l_glo, sst_lp1_glo, diff_sst_glo REAL, DIMENSION(klon_glo):: siv_l_glo, siv_lp1_glo, diff_siv_glo + REAL, DIMENSION(klon_glo,nslay):: bils_glo CHARACTER (len = 20) :: modname = 'limit_slab' + CHARACTER*2 str2 LOGICAL :: read_bils,read_sst,read_siv @@ -44,5 +48,5 @@ ! Initialize saved variables IF (.NOT. ALLOCATED(bils_save)) THEN - ALLOCATE(bils_save(klon), diff_sst_save(klon), diff_siv_save(klon), stat=ierr) + ALLOCATE(bils_save(klon,nslay), diff_sst_save(klon), diff_siv_save(klon), stat=ierr) IF (ierr /= 0) CALL abort_physic('limit_slab', 'pb in allocation',1) END IF @@ -77,11 +81,35 @@ ! ! Read bils_glo - ierr = NF90_INQ_VARID(nid, 'BILS_OCE', nvarid) + bils_glo(:,:)=0. + ! First read first layer + ! try first "BILS_OCE01" + ierr = NF90_INQ_VARID(nid, 'BILS_OCE01', nvarid) IF (ierr /= NF90_NOERR) THEN - read_bils=.FALSE. + ! Else BILS_OCE + ierr = NF90_INQ_VARID(nid, 'BILS_OCE', nvarid) + IF (ierr /= NF90_NOERR) THEN + read_bils=.FALSE. + ELSE + ierr = NF90_GET_VAR(nid,nvarid,bils_glo(:,1),start,epais) + IF (ierr /= NF90_NOERR) read_bils=.FALSE. + ENDIF ELSE - ierr = NF90_GET_VAR(nid,nvarid,bils_glo,start,epais) + ierr = NF90_GET_VAR(nid,nvarid,bils_glo(:,1),start,epais) IF (ierr /= NF90_NOERR) read_bils=.FALSE. END IF + ! Try next layers if more than 1 + IF ((nslay.GT.1).AND.read_bils) THEN + DO i=2,nslay + WRITE(str2,'(i2.2)') i + ierr = NF90_INQ_VARID(nid,'BILS_OCE'//str2, nvarid) + IF (ierr.EQ.NF90_NOERR) THEN + ierr = NF90_GET_VAR(nid,nvarid,bils_glo(:,i),start,epais) + ENDIF + IF (ierr /= NF90_NOERR) THEN + print *,'WARNING : BILS_OCE not found for layer 2' + ENDIF + ENDDO + ENDIF + ! Read sst_glo for this day ierr = NF90_INQ_VARID(nid, 'SST', nvarid) @@ -146,5 +174,5 @@ CALL Scatter(bils_glo, bils_save) ELSE - bils_save(:)=0. + bils_save(:,:)=0. END IF IF (read_sst) THEN @@ -161,5 +189,5 @@ ENDIF ! time to read - lmt_bils(:) = bils_save(:) + lmt_bils(:,:) = bils_save(:,:) diff_sst(:) = diff_sst_save(:) diff_siv(:) = diff_siv_save(:) Index: LMDZ5/trunk/libf/phylmd/ocean_slab_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/ocean_slab_mod.F90 (revision 3001) +++ LMDZ5/trunk/libf/phylmd/ocean_slab_mod.F90 (revision 3002) @@ -1,3 +1,3 @@ -! +!Completed MODULE ocean_slab_mod ! @@ -40,4 +40,10 @@ REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_hdiff !$OMP THREADPRIVATE(dt_hdiff) + ! heat flux convergence due to GM eddy advection + REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_gm + !$OMP THREADPRIVATE(dt_gm) + ! Heat Flux correction + REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_qflux + !$OMP THREADPRIVATE(dt_qflux) ! fraction of ocean covered by sea ice (sic / (oce+sic)) REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: fsic @@ -50,6 +56,4 @@ !$OMP THREADPRIVATE(seaice) ! net surface heat flux, weighted by open ocean fraction - REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: slab_bils - !$OMP THREADPRIVATE(slab_bils) ! slab_bils accumulated over cpl_pas timesteps REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bils_cum @@ -108,4 +112,6 @@ LOGICAL, PUBLIC, SAVE :: slab_hdiff !$OMP THREADPRIVATE(slab_hdiff) + LOGICAL, PUBLIC, SAVE :: slab_gm + !$OMP THREADPRIVATE(slab_gm) REAL, PRIVATE, SAVE :: coef_hdiff ! coefficient for horizontal diffusion !$OMP THREADPRIVATE(coef_hdiff) @@ -155,4 +161,5 @@ ENDIF slabh(1)=50. + CALL getin_p('slab_depth',slabh(1)) IF (nslay.GT.1) THEN slabh(2)=150. @@ -176,4 +183,10 @@ slab_ekman=0 CALL getin_p('slab_ekman',slab_ekman) +! GM eddy advection (2-layers only) + slab_gm=.FALSE. + CALL getin_p('slab_gm',slab_gm) + IF (slab_ekman.LT.2) THEN + slab_gm=.FALSE. + ENDIF ! Convective adjustment IF (nslay.EQ.1) THEN @@ -205,10 +218,4 @@ (modname,'pb allocation tslab', 1) - ALLOCATE(slab_bils(klon), stat = error) - IF (error /= 0) THEN - abort_message='Pb allocation slab_bils' - CALL abort_physic(modname,abort_message,1) - ENDIF - slab_bils(:) = 0.0 ALLOCATE(bils_cum(klon), stat = error) IF (error /= 0) THEN @@ -252,4 +259,20 @@ ENDIF + ALLOCATE(dt_qflux(klon,nslay), stat = error) + IF (error /= 0) THEN + abort_message='Pb allocation dt_qflux' + CALL abort_physic(modname,abort_message,1) + ENDIF + dt_qflux(:,:) = 0.0 + + IF (slab_gm) THEN !GM advection + ALLOCATE(dt_gm(klon,nslay), stat = error) + IF (error /= 0) THEN + abort_message='Pb allocation dt_gm' + CALL abort_physic(modname,abort_message,1) + ENDIF + dt_gm(:,:) = 0.0 + ENDIF + IF (slab_ekman.GT.0) THEN ! ekman transport ALLOCATE(dt_ekman(klon,nslay), stat = error) @@ -276,5 +299,6 @@ IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN CALL gather(zmasq,zmasq_glo) -!$OMP MASTER ! Only master thread +! Master thread/process only +!$OMP MASTER IF (is_mpi_root) THEN CALL ini_slab_transp(zmasq_glo) @@ -317,8 +341,8 @@ radsol, snow, & qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, & - tsurf_new, dflux_s, dflux_l, qflux) + tsurf_new, dflux_s, dflux_l, slab_bils) USE calcul_fluxs_mod - USE slab_heat_transp_mod, ONLY: divgrad_phy, slab_ekman1, slab_ekman2 + USE slab_heat_transp_mod, ONLY: divgrad_phy,slab_ekman1,slab_ekman2,slab_gmdiff USE mod_phys_lmdz_para @@ -369,5 +393,5 @@ REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! new surface tempearture REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l - REAL, DIMENSION(klon), INTENT(OUT) :: qflux + REAL, DIMENSION(klon), INTENT(OUT) :: slab_bils ! Local variables @@ -378,5 +402,6 @@ REAL, DIMENSION(klon) :: cal, beta, dif_grnd ! for Q-Flux computation: d/dt SST, d/dt ice volume (kg/m2), surf fluxes - REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils + REAL, DIMENSION(klon) :: diff_sst, diff_siv + REAL, DIMENSION(klon,nslay) :: lmt_bils ! for surface wind stress REAL, DIMENSION(klon) :: u0, v0 @@ -386,5 +411,5 @@ ! horizontal diffusion and Ekman local vars ! dimension = global domain (klon_glo) instead of // subdomains - REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo, dt_ekman_glo + REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo,dt_ekman_glo,dt_gm_glo ! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_tmp, dt_ekman_tmp @@ -444,9 +469,11 @@ ! !**************************************************************************************** - lmt_bils(:)=0. CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) ! lmt_bils and diff_sst,siv saved by limit_slab - qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400. ! qflux = total QFlux correction (in W/m2) + dt_qflux(:,1)=lmt_bils(:,1)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400. + IF (nslay.GT.1) THEN + dt_qflux(:,2:nslay)=lmt_bils(:,2:nslay) + END IF !**************************************************************************************** @@ -469,6 +496,5 @@ CALL gather(tslab,tslab_glo) ! Compute horiz transport on one process only -!$OMP MASTER ! Only master thread - IF (is_mpi_root) THEN ! Only master processus + IF (is_mpi_root .AND. is_omp_root) THEN ! Only master processus IF (slab_hdiff) THEN dt_hdiff_glo(:,:)=0. @@ -476,4 +502,7 @@ IF (slab_ekman.GT.0) THEN dt_ekman_glo(:,:)=0. + END IF + IF (slab_gm) THEN + dt_gm_glo(:,:)=0. END IF DO i=1,cpl_pas ! time splitting for numerical stability @@ -483,5 +512,5 @@ CALL slab_ekman1(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp) CASE (2) - CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp) + CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp,dt_hdiff_tmp,slab_gm) CASE DEFAULT dt_ekman_tmp(:,:)=0. @@ -493,5 +522,23 @@ ENDDO tslab_glo=tslab_glo+dt_ekman_tmp*dtime + IF (slab_gm) THEN ! Gent-McWilliams eddy advection + dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:) + ! convert dt from K.s-1.(kg.m-2) to K.s-1 + DO k=1,nslay + dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/(slabh(k)*sea_den) + END DO + tslab_glo=tslab_glo+dt_hdiff_tmp*dtime + END IF ENDIF +! GM included in Ekman_2 +! IF (slab_gm) THEN ! Gent-McWilliams eddy advection +! CALL slab_gmdiff(tslab_glo,dt_hdiff_tmp) +! ! convert dt_gm from K.m.s-1 to K.s-1 +! DO k=1,nslay +! dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/slabh(k) +! END DO +! tslab_glo=tslab_glo+dt_hdiff_tmp*dtime +! dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:) +! END IF IF (slab_hdiff) THEN ! horizontal diffusion ! laplacian of slab T @@ -509,10 +556,13 @@ END DO END IF + IF (slab_gm) THEN + !dt_hdiff_glo saved in W/m2 + dt_gm_glo(:,:)=dt_gm_glo(:,:)*sea_cap/cpl_pas + END IF IF (slab_ekman.GT.0) THEN ! dt_ekman_glo saved in W/m2 dt_ekman_glo(:,:)=dt_ekman_glo(:,:)*sea_cap/cpl_pas END IF - END IF ! is_mpi_root -!$OMP END MASTER + END IF ! master process !$OMP BARRIER ! Send new fields back to all processes @@ -520,4 +570,7 @@ IF (slab_hdiff) THEN CALL Scatter(dt_hdiff_glo,dt_hdiff) + END IF + IF (slab_gm) THEN + CALL Scatter(dt_gm_glo,dt_gm) END IF IF (slab_ekman.GT.0) THEN @@ -533,5 +586,8 @@ ! *********************************** ! Add read QFlux - tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas + DO k=1,nslay + tslab(:,k)=tslab(:,k)+dt_qflux(:,k)*cyang*dtime*cpl_pas & + *slabh(1)/slabh(k) + END DO ! Add cumulated surface fluxes tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime @@ -913,5 +969,4 @@ IF (ALLOCATED(tice)) DEALLOCATE(tice) IF (ALLOCATED(seaice)) DEALLOCATE(seaice) - IF (ALLOCATED(slab_bils)) DEALLOCATE(slab_bils) IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg) IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum) @@ -921,4 +976,6 @@ IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman) IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff) + IF (ALLOCATED(dt_gm)) DEALLOCATE(dt_gm) + IF (ALLOCATED(dt_qflux)) DEALLOCATE(dt_qflux) END SUBROUTINE ocean_slab_final Index: LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 3001) +++ LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 3002) @@ -706,4 +706,8 @@ TYPE(ctrl_out), SAVE :: o_tslab = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11/), & 'tslab', 'Temperature ocean slab', 'K', (/ ('', i=1, 10) /)) + TYPE(ctrl_out), SAVE :: o_tslab1 = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11 /), & + 'tslab1', 'Temperature ocean slab', 'K', (/ ('', i=1, 10) /)) + TYPE(ctrl_out), SAVE :: o_tslab2 = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11 /), & + 'tslab2', 'Temperature ocean slab', 'K', (/ ('', i=1, 10) /)) TYPE(ctrl_out), SAVE :: o_slab_tice = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11/), & 'slab_tice', 'Temperature banquise slab', 'K', (/ ('', i=1, 10) /)) @@ -712,4 +716,6 @@ TYPE(ctrl_out), SAVE :: o_slab_hdiff = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11/), & 'slab_hdiff', 'Horizontal diffusion', 'W/m2', (/ ('', i=1, 10) /)) + TYPE(ctrl_out), SAVE :: o_slab_gm = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11 /), & + 'slab_gm', 'GM eddy advection', 'W/m2', (/ ('', i=1, 10) /)) TYPE(ctrl_out), SAVE :: o_slab_ekman = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11/), & 'slab_ekman', 'Ekman heat transport', 'W/m2', (/ ('', i=1, 10) /)) Index: LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90 (revision 3001) +++ LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90 (revision 3002) @@ -91,5 +91,5 @@ o_slab_qflux, o_tslab, o_slab_bils, & o_slab_bilg, o_slab_sic, o_slab_tice, & - o_slab_hdiff, o_slab_ekman, & + o_slab_hdiff, o_slab_ekman, o_slab_gm, & o_weakinv, o_dthmin, o_cldtau, & o_cldemi, o_pr_con_l, o_pr_con_i, & @@ -329,6 +329,6 @@ ok_4xCO2atm - USE ocean_slab_mod, ONLY: nslay, tslab, slab_bils, slab_bilg, tice, & - seaice, slab_ekman,slab_hdiff, dt_ekman, dt_hdiff + USE ocean_slab_mod, ONLY: nslay, tslab, slab_bilg, tice, seaice, & + slab_ekman,slab_hdiff,slab_gm,dt_ekman, dt_hdiff, dt_gm, dt_qflux USE pbl_surface_mod, ONLY: snow USE indice_sol_mod, ONLY: nbsrf @@ -1103,11 +1103,13 @@ ! Output of slab ocean variables IF (type_ocean=='slab ') THEN - CALL histwrite_phy(o_slab_qflux, slab_wfbils) - CALL histwrite_phy(o_slab_bils, slab_bils) + CALL histwrite_phy(o_slab_bils, slab_wfbils) IF (nslay.EQ.1) THEN zx_tmp_fi2d(:)=tslab(:,1) CALL histwrite_phy(o_tslab, zx_tmp_fi2d) + zx_tmp_fi2d(:)=dt_qflux(:,1) + CALL histwrite_phy(o_slab_qflux, zx_tmp_fi2d) ELSE CALL histwrite_phy(o_tslab, tslab(:,1:nslay)) + CALL histwrite_phy(o_slab_qflux, dt_qflux(:,1:nslay)) ENDIF IF (version_ocean=='sicINT') THEN @@ -1116,4 +1118,7 @@ CALL histwrite_phy(o_slab_sic, seaice) ENDIF + IF (slab_gm) THEN + CALL histwrite_phy(o_slab_gm, dt_gm(:,1:nslay)) + END IF IF (slab_hdiff) THEN IF (nslay.EQ.1) THEN Index: LMDZ5/trunk/libf/phylmd/slab_heat_transp_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/slab_heat_transp_mod.F90 (revision 3001) +++ LMDZ5/trunk/libf/phylmd/slab_heat_transp_mod.F90 (revision 3002) @@ -11,4 +11,6 @@ REAL,SAVE,ALLOCATABLE :: fext(:) ! Coriolis f times cell area !$OMP THREADPRIVATE(fext) + REAL,SAVE,ALLOCATABLE :: beta(:) ! df/dy + !$OMP THREADPRIVATE(beta) REAL,SAVE,ALLOCATABLE :: unsairez(:) ! 1/cell area !$OMP THREADPRIVATE(unsairez) @@ -34,10 +36,26 @@ !$OMP THREADPRIVATE(airev) - ! Local variables for horiz mass flux in slab - LOGICAL,SAVE :: alpha_var + ! Local parameters for slab transport + LOGICAL,SAVE :: alpha_var ! variable coef for deep temp (1 layer) !$OMP THREADPRIVATE(alpha_var) - LOGICAL,SAVE :: slab_upstream + LOGICAL,SAVE :: slab_upstream ! upstream scheme ? (1 layer) !$OMP THREADPRIVATE(slab_upstream) - + LOGICAL,SAVE :: slab_sverdrup ! use wind stress curl at equator + !$OMP THREADPRIVATE(slab_sverdrup) + LOGICAL,SAVE :: slab_tyeq ! use merid wind stress at equator + !$OMP THREADPRIVATE(slab_tyeq) + LOGICAL,SAVE :: ekman_zonadv ! use zonal advection by Ekman currents + !$OMP THREADPRIVATE(ekman_zonadv) + LOGICAL,SAVE :: ekman_zonavg ! zonally average wind stress + !$OMP THREADPRIVATE(ekman_zonavg) + + REAL,SAVE :: alpham + !$OMP THREADPRIVATE(alpham) + REAL,SAVE :: gmkappa + !$OMP THREADPRIVATE(gmkappa) + REAL,SAVE :: gm_smax + !$OMP THREADPRIVATE(gm_smax) + +! geometry variables : f, beta, mask... REAL,SAVE,ALLOCATABLE :: zmasqu(:) ! continent mask for zonal mass flux !$OMP THREADPRIVATE(zmasqu) @@ -46,4 +64,6 @@ REAL,SAVE,ALLOCATABLE :: unsfv(:) ! 1/f, v points !$OMP THREADPRIVATE(unsfv) + REAL,SAVE,ALLOCATABLE :: unsbv(:) ! 1/beta + !$OMP THREADPRIVATE(unsbv) REAL,SAVE,ALLOCATABLE :: unsev(:) ! 1/epsilon (drag) !$OMP THREADPRIVATE(unsev) @@ -63,5 +83,6 @@ cu_,cuvsurcv_,cv_,cvusurcu_, & aire_,apoln_,apols_, & - aireu_,airev_) + aireu_,airev_,rlatv) + USE comconst_mod, ONLY: omeg, rad ! number of points in lon, lat IMPLICIT NONE @@ -82,5 +103,5 @@ REAL,INTENT(IN) :: aireu_(ip1jmp1) REAL,INTENT(IN) :: airev_(ip1jm) - + REAL,INTENT(IN) :: rlatv(nbp_lat-1) ! Sanity check on dimensions @@ -112,5 +133,7 @@ aireu(:)=aireu_(:) allocate(airev(ip1jm)) - airev(:)=airev_(:) + airev(:)=airev_(:) + allocate(beta(nbp_lat-1)) + beta(:)=2*omeg*cos(rlatv(:))/rad END SUBROUTINE ini_slab_transp_geom @@ -139,4 +162,42 @@ ip1jmp1=(nbp_lon+1)*nbp_lat +! Options for Heat transport + ! Alpha variable? + alpha_var=.FALSE. + CALL getin('slab_alphav',alpha_var) + print *,'alpha variable',alpha_var +! centered ou upstream scheme for meridional transport + slab_upstream=.FALSE. + CALL getin('slab_upstream',slab_upstream) + print *,'upstream slab scheme',slab_upstream +! Sverdrup balance at equator ? + slab_sverdrup=.TRUE. + CALL getin('slab_sverdrup',slab_sverdrup) + print *,'Sverdrup balance',slab_sverdrup +! Use tauy for meridional flux at equator ? + slab_tyeq=.TRUE. + CALL getin('slab_tyeq',slab_tyeq) + print *,'Tauy forcing at equator',slab_tyeq +! Use tauy for meridional flux at equator ? + ekman_zonadv=.TRUE. + CALL getin('slab_ekman_zonadv',ekman_zonadv) + print *,'Use Ekman flow in zonal direction',ekman_zonadv +! Use tauy for meridional flux at equator ? + ekman_zonavg=.FALSE. + CALL getin('slab_ekman_zonavg',ekman_zonavg) + print *,'Use zonally-averaged wind stress ?',ekman_zonavg +! Value of alpha + alpham=2./3. + CALL getin('slab_alpha',alpham) + print *,'slab_alpha',alpham +! GM k coefficient (m2/s) for 2-layers + gmkappa=1000. + CALL getin('slab_gmkappa',gmkappa) + print *,'slab_gmkappa',gmkappa +! GM k coefficient (m2/s) for 2-layers + gm_smax=2e-3 + CALL getin('slab_gm_smax',gm_smax) + print *,'slab_gm_smax',gm_smax +! ----------------------------------------------------------- ! Define ocean / continent mask (no flux into continent cell) allocate(zmasqu(ip1jmp1)) @@ -161,7 +222,9 @@ END IF END DO + +! ----------------------------------------------------------- +! Coriolis and friction for Ekman transport slab_ekman=2 CALL getin("slab_ekman",slab_ekman) -! Coriolis and friction for Ekman transport IF (slab_ekman.GT.0) THEN allocate(unsfv(ip1jm)) @@ -169,4 +232,5 @@ allocate(unsfu(ip1jmp1)) allocate(unseu(ip1jmp1)) + allocate(unsbv(ip1jm)) eps=1e-5 ! Drag @@ -176,22 +240,30 @@ ! coefs to convert tau_x, tau_y to Ekman mass fluxes ! on 2D grid v points + ! Compute correction factor [0 1] near the equator (f<