MODULE dynadv_tam #ifdef key_tam !!============================================================================== !! *** MODULE dynadv_tam *** !! Ocean active tracers: advection scheme control !! Tangent and Adjoint module !!============================================================================== !! History of the direct module: !! 9.0 ! 06-11 (G. Madec) Original code !! History of the TAM module: !! 9.0 ! 08-08 (A. Vidard) first version !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! dyn_adv : compute the momentum advection trend !! dyn_adv_ctl : control the different options of advection scheme !!---------------------------------------------------------------------- USE par_kind, ONLY: & ! Precision variables & wp USE par_oce, ONLY: & ! Ocean space and time domain variables & jpi, & & jpj, & & jpk, & & jpiglo USE oce, ONLY: & & un, & & vn, & & wn, & & ua, & & va USE dom_oce, ONLY: & ! ocean space and time domain & umask, & & vmask, & & mig, & & mjg, & & e1u, & & e2u, & & e1v, & & e2v, & # if defined key_zco & e3t_0, & # else & e3u, & & e3v, & # endif & nldi, & & nldj, & & nlei, & & nlej USE oce_tam, ONLY: & & un_tl, & & vn_tl, & & wn_tl, & & ua_tl, & & va_tl, & & un_ad, & & vn_ad, & & wn_ad, & & ua_ad, & & va_ad USE in_out_manager, ONLY: & ! I/O manager & ctl_stop, & & lwp, & & lk_esopa, & & numout, & & numnam, & & nit000, & & nitend USE gridrandom, ONLY: & ! Random Gaussian noise on grids & grid_random USE dotprodfld, ONLY: & ! Computes dot product for 3D and 2D fields & dot_product USE tstool_tam, ONLY: & & prntst_adj, & ! & prntst_tlm, & ! & stdu, & ! stdev for u-velocity & stdv, & ! v-velocity & stdw ! w-velocity USE dynadv, ONLY: & & dyn_adv_ctl, & & ln_dynadv_vec, & ! vector form flag & ln_dynadv_cen2, & ! flux form - 2nd order centered scheme flag & ln_dynadv_ubs ! flux form - 3rd order UBS s ! USE dynadv_cen2_tam ! centred flux form advection (dyn_adv_cen2 routine) ! USE dynadv_ubs_tam ! UBS flux form advection (dyn_adv_ubs routine) USE dynkeg_tam, ONLY: & ! kinetic energy gradient (dyn_keg routine) & dyn_keg_tan, & & dyn_keg_adj USE dynzad_tam, ONLY: & ! vertical advection (dyn_zad routine) & dyn_zad_tan, & & dyn_zad_adj USE wzvmod , ONLY: & & wzv USE divcur , ONLY: & & div_cur USE wzvmod_tam , ONLY: & & wzv_tan USE divcur_tam , ONLY: & & div_cur_tan IMPLICIT NONE PRIVATE PUBLIC dyn_adv_tan ! routine called by steptan module PUBLIC dyn_adv_adj ! routine called by stepadj module PUBLIC dyn_adv_adj_tst ! routine called by the tst module PUBLIC dyn_adv_tlm_tst ! routine called by tamtst.F90 INTEGER :: nadv ! choice of the formulation and scheme for the advection LOGICAL :: lfirst=.TRUE. !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" CONTAINS SUBROUTINE dyn_adv_tan( kt ) !!--------------------------------------------------------------------- !! *** ROUTINE dyn_adv_tan *** !! !! ** Purpose of the direct routine: !! compute the ocean momentum advection trend. !! !! ** Method : - Update (ua,va) with the advection term following nadv !! NB: in flux form advection (ln_dynadv_cen2 or ln_dynadv_ubs=T) !! a metric term is add to the coriolis term while in vector form !! it is the relative vorticity which is added to coriolis term !! (see dynvor module). !!---------------------------------------------------------------------- INTEGER, INTENT( in ) :: kt ! ocean time-step index !!---------------------------------------------------------------------- ! IF( kt == nit000 ) CALL dyn_adv_ctl_tam! initialisation & control of options SELECT CASE ( nadv ) ! compute advection trend and add it to general trend CASE ( 0 ) CALL dyn_keg_tan ( kt ) ! vector form : horizontal gradient of kinetic energy CALL dyn_zad_tan ( kt ) ! vector form : vertical advection CASE ( 1 ) IF (lwp) WRITE(numout,*) 'dyn_adv_cen2_tan not available yet' CALL abort ! CALL dyn_adv_cen2_tan( kt ) ! 2nd order centered scheme CASE ( 2 ) IF (lwp) WRITE(numout,*) 'dyn_adv_ubs_tan not available yet' CALL abort ! CALL dyn_adv_ubs_tan ( kt ) ! 3rd order UBS scheme ! CASE (-1 ) ! esopa: test all possibility with control print CALL dyn_keg_tan ( kt ) CALL dyn_zad_tan ( kt ) ! CALL dyn_adv_cen2_tan( kt ) ! CALL dyn_adv_ubs_tan ( kt ) END SELECT ! END SUBROUTINE dyn_adv_tan SUBROUTINE dyn_adv_adj( kt ) !!--------------------------------------------------------------------- !! *** ROUTINE dyn_adv_adj *** !! !! ** Purpose of the direct routine: !! compute the ocean momentum advection trend. !! !! ** Method : - Update (ua,va) with the advection term following nadv !! NB: in flux form advection (ln_dynadv_cen2 or ln_dynadv_ubs=T) !! a metric term is add to the coriolis term while in vector form !! it is the relative vorticity which is added to coriolis term !! (see dynvor module). !!---------------------------------------------------------------------- INTEGER, INTENT( in ) :: kt ! ocean time-step index !!---------------------------------------------------------------------- ! IF ( kt == nitend ) CALL dyn_adv_ctl_tam SELECT CASE ( nadv ) ! compute advection trend and add it to general trend CASE ( 0 ) CALL dyn_zad_adj ( kt ) ! vector form : vertical advection CALL dyn_keg_adj ( kt ) ! vector form : horizontal gradient of kinetic energy CASE ( 1 ) IF (lwp) WRITE(numout,*) 'dyn_adv_cen2_adj not available yet' CALL abort ! CALL dyn_adv_cen2_adj( kt ) ! 2nd order centered scheme CASE ( 2 ) IF (lwp) WRITE(numout,*) 'dyn_adv_ubs_adj not available yet' CALL abort ! CALL dyn_adv_ubs_adj ( kt ) ! 3rd order UBS scheme ! CASE (-1 ) ! esopa: test all possibility with control print ! CALL dyn_adv_ubs_adj ( kt ) ! CALL dyn_adv_cen2_adj( kt ) CALL dyn_zad_adj ( kt ) CALL dyn_keg_adj ( kt ) END SELECT ! END SUBROUTINE dyn_adv_adj SUBROUTINE dyn_adv_adj_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE dyn_adv_adj_tst *** !! !! ** Purpose : Test the adjoint routine. !! !! ** Method : Verify the scalar product !! !! ( L dx )^T W dy = dx^T L^T W dy !! !! where L = tangent routine !! L^T = adjoint routine !! W = diagonal matrix of scale factors !! dx = input perturbation (random field) !! dy = L dx !! !! !! History : !! ! 08-08 (A. Vidard) !!----------------------------------------------------------------------- !! * Modules used !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit !! * Local declarations INTEGER :: & & ji, & ! dummy loop indices & jj, & & jk INTEGER, DIMENSION(jpi,jpj) :: & & iseed_2d ! 2D seed for the random number generator REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zun_tlin, & ! Tangent input: now u-velocity & zvn_tlin, & ! Tangent input: now v-velocity & zwn_tlin, & ! Tangent input: now w-velocity & zua_tlin, & ! Tangent input: after u-velocity & zva_tlin, & ! Tangent input: after u-velocity & zua_tlout, & ! Tangent output:after u-velocity & zva_tlout, & ! Tangent output:after v-velocity & zua_adin, & ! adjoint input: after u-velocity & zva_adin, & ! adjoint input: after v-velocity & zun_adout, & ! adjoint output: now u-velocity & zvn_adout, & ! adjoint output: now v-velocity & zwn_adout, & ! adjoint output: now u-velocity & zua_adout, & ! adjoint output:after v-velocity & zva_adout, & ! adjoint output:after u-velocity & zuvw ! 3D random field for u, v and w REAL(KIND=wp) :: & & zsp1, & ! scalar product involving the tangent routine & zsp1_1, & ! scalar product components & zsp1_2, & & zsp2, & ! scalar product involving the adjoint routine & zsp2_1, & ! scalar product components & zsp2_2, & & zsp2_3, & & zsp2_4, & & zsp2_5 CHARACTER(LEN=14) :: cl_name ! Allocate memory ALLOCATE( & & zun_tlin(jpi,jpj,jpk), & & zvn_tlin(jpi,jpj,jpk), & & zwn_tlin(jpi,jpj,jpk), & & zua_tlin(jpi,jpj,jpk), & & zva_tlin(jpi,jpj,jpk), & & zua_tlout(jpi,jpj,jpk), & & zva_tlout(jpi,jpj,jpk), & & zua_adin(jpi,jpj,jpk), & & zva_adin(jpi,jpj,jpk), & & zun_adout(jpi,jpj,jpk), & & zvn_adout(jpi,jpj,jpk), & & zwn_adout(jpi,jpj,jpk), & & zua_adout(jpi,jpj,jpk), & & zva_adout(jpi,jpj,jpk), & & zuvw(jpi,jpj,jpk) & & ) !================================================================== ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and ! dy = ( hdivb_tl, hdivn_tl ) !================================================================== !-------------------------------------------------------------------- ! Reset the tangent and adjoint variables !-------------------------------------------------------------------- zun_tlin(:,:,:) = 0.0_wp zvn_tlin(:,:,:) = 0.0_wp zwn_tlin(:,:,:) = 0.0_wp zua_tlin(:,:,:) = 0.0_wp zva_tlin(:,:,:) = 0.0_wp zua_tlout(:,:,:) = 0.0_wp zva_tlout(:,:,:) = 0.0_wp zua_adin(:,:,:) = 0.0_wp zva_adin(:,:,:) = 0.0_wp zun_adout(:,:,:) = 0.0_wp zvn_adout(:,:,:) = 0.0_wp zwn_adout(:,:,:) = 0.0_wp zua_adout(:,:,:) = 0.0_wp zva_adout(:,:,:) = 0.0_wp zuvw(:,:,:) = 0.0_wp un_tl(:,:,:) = 0.0_wp vn_tl(:,:,:) = 0.0_wp wn_tl(:,:,:) = 0.0_wp ua_tl(:,:,:) = 0.0_wp va_tl(:,:,:) = 0.0_wp un_ad(:,:,:) = 0.0_wp vn_ad(:,:,:) = 0.0_wp wn_ad(:,:,:) = 0.0_wp ua_ad(:,:,:) = 0.0_wp va_ad(:,:,:) = 0.0_wp !recompute wn from un and vn CALL div_cur( nit000 ) CALL wzv( nit000 ) !-------------------------------------------------------------------- ! Initialize the tangent input with random noise: dx !-------------------------------------------------------------------- DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 596035 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, zuvw, 'U', 0.0_wp, stdu ) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zun_tlin(ji,jj,jk) = zuvw(ji,jj,jk) END DO END DO END DO DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 523432 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, zuvw, 'V', 0.0_wp, stdv ) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zvn_tlin(ji,jj,jk) = zuvw(ji,jj,jk) END DO END DO END DO DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 456953 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, zuvw, 'W', 0.0_wp, stdw ) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zwn_tlin(ji,jj,jk) = zuvw(ji,jj,jk) END DO END DO END DO DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 432545 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, zuvw, 'U', 0.0_wp, stdu ) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zua_tlin(ji,jj,jk) = zuvw(ji,jj,jk) END DO END DO END DO DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 287503 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, zuvw, 'V', 0.0_wp, stdv ) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zva_tlin(ji,jj,jk) = zuvw(ji,jj,jk) END DO END DO END DO un_tl(:,:,:) = zun_tlin(:,:,:) vn_tl(:,:,:) = zvn_tlin(:,:,:) !recompute wn_tl from un_tl and vn_tl CALL div_cur_tan( nit000 ) CALL wzv_tan( nit000 ) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zwn_tlin(ji,jj,jk) = wn_tl(ji,jj,jk) END DO END DO END DO wn_tl(:,:,:) = zwn_tlin(:,:,:) ua_tl(:,:,:) = zua_tlin(:,:,:) va_tl(:,:,:) = zva_tlin(:,:,:) CALL dyn_adv_tan ( nit000 ) zua_tlout(:,:,:) = ua_tl(:,:,:) zva_tlout(:,:,:) = va_tl(:,:,:) !-------------------------------------------------------------------- ! Initialize the adjoint variables: dy^* = W dy !-------------------------------------------------------------------- DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zua_adin(ji,jj,jk) = zua_tlout(ji,jj,jk) & & * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) & & * umask(ji,jj,jk) zva_adin(ji,jj,jk) = zva_tlout(ji,jj,jk) & & * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) & & * vmask(ji,jj,jk) END DO END DO END DO !-------------------------------------------------------------------- ! Compute the scalar product: ( L dx )^T W dy !-------------------------------------------------------------------- zsp1_1 = DOT_PRODUCT( zua_tlout, zua_adin ) zsp1_2 = DOT_PRODUCT( zva_tlout, zva_adin ) zsp1 = zsp1_1 + zsp1_2 !-------------------------------------------------------------------- ! Call the adjoint routine: dx^* = L^T dy^* !-------------------------------------------------------------------- ua_ad(:,:,:) = zua_adin(:,:,:) va_ad(:,:,:) = zva_adin(:,:,:) CALL dyn_adv_adj ( nit000 ) zun_adout(:,:,:) = un_ad(:,:,:) zvn_adout(:,:,:) = vn_ad(:,:,:) zwn_adout(:,:,:) = wn_ad(:,:,:) zua_adout(:,:,:) = ua_ad(:,:,:) zva_adout(:,:,:) = va_ad(:,:,:) zsp2_1 = DOT_PRODUCT( zun_tlin, zun_adout ) zsp2_2 = DOT_PRODUCT( zvn_tlin, zvn_adout ) zsp2_3 = DOT_PRODUCT( zwn_tlin, zwn_adout ) zsp2_4 = DOT_PRODUCT( zua_tlin, zua_adout ) zsp2_5 = DOT_PRODUCT( zva_tlin, zva_adout ) zsp2 = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4 + zsp2_5 ! Compare the scalar products cl_name = 'dyn_adv_adj ' CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 ) DEALLOCATE( & & zun_tlin, & & zvn_tlin, & & zwn_tlin, & & zua_tlin, & & zva_tlin, & & zua_tlout, & & zva_tlout, & & zua_adin, & & zva_adin, & & zun_adout, & & zvn_adout, & & zwn_adout, & & zua_adout, & & zva_adout, & & zuvw & & ) END SUBROUTINE dyn_adv_adj_tst !!====================================================================== SUBROUTINE dyn_adv_ctl_tam !!--------------------------------------------------------------------- !! *** ROUTINE dyn_adv_ctl *** !! !! ** Purpose : Control the consistency between namelist options for !! momentum advection formulation & scheme and set nadv !!---------------------------------------------------------------------- INTEGER :: ioptio NAMELIST/nam_dynadv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs !!---------------------------------------------------------------------- IF (lfirst) THEN REWIND ( numnam ) ! Read Namelist nam_dynadv : momentum advection scheme READ ( numnam, nam_dynadv ) IF(lwp) THEN ! Namelist print WRITE(numout,*) WRITE(numout,*) 'dyn_adv_ctl : choice/control of the momentum advection scheme' WRITE(numout,*) '~~~~~~~~~~~' WRITE(numout,*) ' Namelist nam_dynadv : chose a advection formulation & scheme for momentum' WRITE(numout,*) ' Vector/flux form (T/F) ln_dynadv_vec = ', ln_dynadv_vec WRITE(numout,*) ' 2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2 WRITE(numout,*) ' 3rd order UBS advection scheme ln_dynadv_ubs = ', ln_dynadv_ubs ENDIF ioptio = 0 ! Parameter control IF( ln_dynadv_vec ) ioptio = ioptio + 1 IF( ln_dynadv_cen2 ) ioptio = ioptio + 1 IF( ln_dynadv_ubs ) ioptio = ioptio + 1 IF( lk_esopa ) ioptio = 1 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE advection scheme in namelist nam_dynadv' ) ! ! Set nadv IF( ln_dynadv_vec ) nadv = 0 IF( ln_dynadv_cen2 ) nadv = 1 IF( ln_dynadv_ubs ) nadv = 2 IF( lk_esopa ) nadv = -1 IF(lwp) THEN ! Print the choice WRITE(numout,*) IF( nadv == 0 ) WRITE(numout,*) ' vector form : keg + zad + vor is used' IF( nadv == 1 ) WRITE(numout,*) ' flux form : 2nd order scheme is used' IF( nadv == 2 ) WRITE(numout,*) ' flux form : UBS scheme is used' IF( nadv == -1 ) WRITE(numout,*) ' esopa test: use all advection formulation' ENDIF ! lfirst = .FALSE. END IF END SUBROUTINE dyn_adv_ctl_tam #if defined key_tst_tlm SUBROUTINE dyn_adv_tlm_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE dyn_adv_tlm_tst *** !! !! ** Purpose : Test the adjoint routine. !! !! ** Method : Verify the tangent with Taylor expansion !! !! M(x+hdx) = M(x) + L(hdx) + O(h^2) !! !! where L = tangent routine !! M = direct routine !! dx = input perturbation (random field) !! h = ration on perturbation !! !! In the tangent test we verify that: !! M(x+h*dx) - M(x) !! g(h) = ------------------ ---> 1 as h ---> 0 !! L(h*dx) !! and !! g(h) - 1 !! f(h) = ---------- ---> k (costant) as h ---> 0 !! p !! !! History : !! ! 09-08 (A. Vigilant) !!----------------------------------------------------------------------- !! * Modules used USE dynadv USE tamtrj ! writing out state trajectory USE par_tlm, ONLY: & & tlm_bch, & & cur_loop, & & h_ratio USE istate_mod USE divcur ! horizontal divergence and relative vorticity USE wzvmod ! vertical velocity USE gridrandom, ONLY: & & grid_rd_sd USE trj_tam USE oce , ONLY: & ! ocean dynamics and tracers variables & ua, va USE opatam_tst_ini, ONLY: & & tlm_namrd USE tamctl, ONLY: & ! Control parameters & numtan, numtan_sc !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit !! * Local declarations INTEGER :: & & ji, & ! dummy loop indices & jj, & & jk REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zun_tlin, & ! Tangent input: now u-velocity & zvn_tlin, & ! Tangent input: now v-velocity & zwn_tlin, & ! Tangent input: now w-velocity & zua_tlin, & ! Tangent input: after u-velocity & zva_tlin, & ! Tangent input: after u-velocity & zua_out, & ! Tangent output:after u-velocity & zva_out, & ! Tangent output:after v-velocity & zua_wop, & ! Tangent output:after u-velocity & zva_wop, & ! Tangent output:after v-velocity & z3r ! 3D field REAL(KIND=wp) :: & & zsp1, & ! scalar product & zsp1_1, & ! scalar product & zsp1_2, & & zsp2, & ! scalar product & zsp2_1, & ! scalar product & zsp2_2, & & zsp3, & ! scalar product & zsp3_1, & ! scalar product & zsp3_2, & & zsp2_3, & & zsp2_4, & & zzsp, & & zzsp_1, & & zzsp_2, & & gamma, & & zgsp1, & & zgsp2, & & zgsp3, & & zgsp4, & & zgsp5, & & zgsp6, & & zgsp7 CHARACTER(LEN=14) :: cl_name CHARACTER (LEN=128) :: file_out, file_wop, file_xdx CHARACTER (LEN=90) :: FMT REAL(KIND=wp), DIMENSION(100):: & & zscua, zscva, & & zscerrua, & & zscerrva INTEGER, DIMENSION(100):: & & iiposua, iiposva, & & ijposua, ijposva, & & ikposua, ikposva INTEGER:: & & ii, & & isamp=40, & & jsamp=40, & & ksamp=10, & & numsctlm REAL(KIND=wp), DIMENSION(jpi,jpj,jpk) :: & & zerrua, zerrva ! Allocate memory ALLOCATE( & & zun_tlin(jpi,jpj,jpk), & & zvn_tlin(jpi,jpj,jpk), & & zwn_tlin(jpi,jpj,jpk), & & zua_tlin(jpi,jpj,jpk), & & zva_tlin(jpi,jpj,jpk), & & zua_out(jpi,jpj,jpk), & & zva_out(jpi,jpj,jpk), & & zua_wop(jpi,jpj,jpk), & & zva_wop(jpi,jpj,jpk), & & z3r (jpi,jpj,jpk) & & ) !-------------------------------------------------------------------- ! Reset variables !-------------------------------------------------------------------- zun_tlin( :,:,:) = 0.0_wp zvn_tlin( :,:,:) = 0.0_wp zwn_tlin( :,:,:) = 0.0_wp zua_tlin( :,:,:) = 0.0_wp zva_tlin( :,:,:) = 0.0_wp zua_out ( :,:,:) = 0.0_wp zva_out ( :,:,:) = 0.0_wp zua_wop ( :,:,:) = 0.0_wp zva_wop ( :,:,:) = 0.0_wp zscua(:) = 0.0_wp zscva(:) = 0.0_wp zscerrua(:) = 0.0_wp zscerrva(:) = 0.0_wp zerrua(:,:,:) = 0.0_wp zerrva(:,:,:) = 0.0_wp !-------------------------------------------------------------------- ! Output filename Xn=F(X0) !-------------------------------------------------------------------- CALL tlm_namrd gamma = h_ratio file_wop='trj_wop_dynadv' file_xdx='trj_xdx_dynadv' !-------------------------------------------------------------------- ! Initialize the tangent input with random noise: dx !-------------------------------------------------------------------- IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN CALL grid_rd_sd( 596035, z3r, 'U', 0.0_wp, stdu) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zun_tlin(ji,jj,jk) = z3r(ji,jj,jk) END DO END DO END DO CALL grid_rd_sd( 523432, z3r, 'V', 0.0_wp, stdv) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zvn_tlin(ji,jj,jk) = z3r(ji,jj,jk) END DO END DO END DO CALL grid_rd_sd( 432545, z3r, 'U', 0.0_wp, stdu) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zua_tlin(ji,jj,jk) = z3r(ji,jj,jk) END DO END DO END DO CALL grid_rd_sd( 287503, z3r, 'V', 0.0_wp, stdv) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zva_tlin(ji,jj,jk) = z3r(ji,jj,jk) END DO END DO END DO ENDIF !-------------------------------------------------------------------- ! Complete Init for Direct !------------------------------------------------------------------- IF ( tlm_bch /= 2 ) CALL istate_p ! *** initialize the reference trajectory ! ------------ CALL trj_rea( nit000-1, 1 ) CALL trj_rea( nit000, 1 ) ua(:,:,:) = un(:,:,:) va(:,:,:) = vn(:,:,:) IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN zun_tlin(:,:,:) = gamma * zun_tlin(:,:,:) un(:,:,:) = un(:,:,:) + zun_tlin(:,:,:) zvn_tlin(:,:,:) = gamma * zvn_tlin(:,:,:) vn(:,:,:) = vn(:,:,:) + zvn_tlin(:,:,:) zua_tlin(:,:,:) = gamma * zua_tlin(:,:,:) ua(:,:,:) = ua(:,:,:) + zua_tlin(:,:,:) zva_tlin(:,:,:) = gamma * zva_tlin(:,:,:) va(:,:,:) = va(:,:,:) + zva_tlin(:,:,:) !recompute wn from un and vn CALL div_cur( nit000) CALL wzv( nit000) ENDIF !-------------------------------------------------------------------- ! Compute the direct model F(X0,t=n) = Xn !-------------------------------------------------------------------- IF ( tlm_bch /= 2 ) CALL dyn_adv(nit000) IF ( tlm_bch == 0 ) CALL trj_wri_spl(file_wop) IF ( tlm_bch == 1 ) CALL trj_wri_spl(file_xdx) !-------------------------------------------------------------------- ! Compute the Tangent !-------------------------------------------------------------------- IF ( tlm_bch == 2 ) THEN !-------------------------------------------------------------------- ! Initialize the tangent variables !-------------------------------------------------------------------- CALL trj_rea( nit000-1, 1 ) CALL trj_rea( nit000, 1 ) ua(:,:,:) = un(:,:,:) va(:,:,:) = vn(:,:,:) un_tl (:,:,:) = zun_tlin (:,:,:) vn_tl (:,:,:) = zvn_tlin (:,:,:) !recompute wn_tl from un_tl and vn_tl CALL div_cur_tan( nit000 ) CALL wzv_tan( nit000 ) ua_tl (:,:,:) = zua_tlin (:,:,:) va_tl (:,:,:) = zva_tlin (:,:,:) CALL dyn_adv_tan(nit000) !-------------------------------------------------------------------- ! Compute the scalar product: ( L(t0,tn) gamma dx0 ) ) !-------------------------------------------------------------------- zsp2_1 = DOT_PRODUCT( ua_tl, ua_tl ) zsp2_2 = DOT_PRODUCT( va_tl, va_tl ) zsp2 = zsp2_1 + zsp2_2 !-------------------------------------------------------------------- ! Storing data !-------------------------------------------------------------------- CALL trj_rd_spl(file_wop) zua_wop (:,:,:) = ua (:,:,:) zva_wop (:,:,:) = va (:,:,:) CALL trj_rd_spl(file_xdx) zua_out (:,:,:) = ua (:,:,:) zva_out (:,:,:) = va (:,:,:) !-------------------------------------------------------------------- ! Compute the Linearization Error ! Nn = M( X0+gamma.dX0, t0,tn) - M(X0, t0,tn) ! and ! Compute the Linearization Error ! En = Nn -TL(gamma.dX0, t0,tn) !-------------------------------------------------------------------- ! Warning: Here we re-use local variables z()_out and z()_wop ii=0 DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi zua_out (ji,jj,jk) = zua_out (ji,jj,jk) - zua_wop (ji,jj,jk) zua_wop (ji,jj,jk) = zua_out (ji,jj,jk) - ua_tl (ji,jj,jk) IF ( ua_tl(ji,jj,jk) .NE. 0.0_wp ) & & zerrua(ji,jj,jk) = zua_out(ji,jj,jk)/ua_tl(ji,jj,jk) IF( (MOD(ji, isamp) .EQ. 0) .AND. & & (MOD(jj, jsamp) .EQ. 0) .AND. & & (MOD(jk, ksamp) .EQ. 0) ) THEN ii = ii+1 iiposua(ii) = ji ijposua(ii) = jj ikposua(ii) = jk IF ( INT(umask(ji,jj,jk)) .NE. 0) THEN zscua (ii) = zua_wop(ji,jj,jk) zscerrua (ii) = ( zerrua(ji,jj,jk) - 1.0_wp ) / gamma ENDIF ENDIF END DO END DO END DO ii=0 DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi zva_out (ji,jj,jk) = zva_out (ji,jj,jk) - zva_wop (ji,jj,jk) zva_wop (ji,jj,jk) = zva_out (ji,jj,jk) - va_tl (ji,jj,jk) IF ( va_tl(ji,jj,jk) .NE. 0.0_wp ) & & zerrva(ji,jj,jk) = zva_out(ji,jj,jk)/va_tl(ji,jj,jk) IF( (MOD(ji, isamp) .EQ. 0) .AND. & & (MOD(jj, jsamp) .EQ. 0) .AND. & & (MOD(jk, ksamp) .EQ. 0) ) THEN ii = ii+1 iiposva(ii) = ji ijposva(ii) = jj ikposva(ii) = jk IF ( INT(vmask(ji,jj,jk)) .NE. 0) THEN zscva (ii) = zua_wop(ji,jj,jk) zscerrva (ii) = ( zerrva(ji,jj,jk) - 1.0_wp ) / gamma ENDIF ENDIF END DO END DO END DO zsp1_1 = DOT_PRODUCT( zua_out, zua_out ) zsp1_2 = DOT_PRODUCT( zva_out, zva_out ) zsp1 = zsp1_1 + zsp1_2 zsp3_1 = DOT_PRODUCT( zua_wop, zua_wop ) zsp3_2 = DOT_PRODUCT( zva_wop, zva_wop ) zsp3 = zsp3_1 + zsp3_2 !-------------------------------------------------------------------- ! Print the linearization error En - norme 2 !-------------------------------------------------------------------- ! 14 char:'12345678901234' cl_name = 'dynadv_tam:En ' zzsp = SQRT(zsp3) zzsp_1 = SQRT(zsp3_1) zzsp_2 = SQRT(zsp3_2) zgsp5 = zzsp CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) !-------------------------------------------------------------------- ! Compute TLM norm2 !-------------------------------------------------------------------- zzsp = SQRT(zsp2) zzsp_1 = SQRT(zsp2_1) zzsp_2 = SQRT(zsp2_2) zgsp4 = zzsp cl_name = 'dynadv_tam:Ln2' CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) !-------------------------------------------------------------------- ! Print the linearization error Nn - norme 2 !-------------------------------------------------------------------- zzsp = SQRT(zsp1) zzsp_1 = SQRT(zsp1_1) zzsp_2 = SQRT(zsp1_2) cl_name = 'dynadv:Mhdx-Mx' CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) zgsp3 = SQRT( zsp3/zsp2 ) zgsp7 = zgsp3/gamma zgsp1 = zzsp zgsp2 = zgsp1 / zgsp4 zgsp6 = (zgsp2 - 1.0_wp)/gamma FMT = "(A8,2X,I4.4,2X,E6.1,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13)" WRITE(numtan,FMT) 'dynadv ', cur_loop, h_ratio, zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, zgsp6, zgsp7 !-------------------------------------------------------------------- ! Unitary calculus !-------------------------------------------------------------------- FMT = "(A8,2X,A8,2X,I4.4,2X,E6.1,2X,I4.4,2X,I4.4,2X,I4.4,2X,E20.13,1X)" cl_name = 'dynadv ' IF (lwp) THEN DO ii=1, 100, 1 IF ( zscua(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscua ', & & cur_loop, h_ratio, ii, iiposua(ii), ijposua(ii), zscua(ii) ENDDO DO ii=1, 100, 1 IF ( zscva(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscva ', & & cur_loop, h_ratio, ii, iiposva(ii), ijposva(ii), zscva(ii) ENDDO DO ii=1, 100, 1 IF ( zscerrua(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrua ', & & cur_loop, h_ratio, ii, iiposua(ii), ijposua(ii), zscerrua(ii) ENDDO DO ii=1, 100, 1 IF ( zscerrva(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrva ', & & cur_loop, h_ratio, ii, iiposva(ii), ijposva(ii), zscerrva(ii) ENDDO ! write separator WRITE(numtan_sc,"(A4)") '====' ENDIF ENDIF DEALLOCATE( & & zun_tlin, zvn_tlin, zwn_tlin, & & zua_tlin, zva_tlin, & & zua_out, zva_out, & & zua_wop, zva_wop, & & z3r & & ) END SUBROUTINE dyn_adv_tlm_tst #endif #endif END MODULE dynadv_tam