Changeset 10969 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ASM/asminc.F90
- Timestamp:
- 2019-05-13T12:02:13+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ASM/asminc.F90
r10954 r10969 102 102 CONTAINS 103 103 104 SUBROUTINE asm_inc_init( K mm)104 SUBROUTINE asm_inc_init( Kbb, Kmm, Krhs ) 105 105 !!---------------------------------------------------------------------- 106 106 !! *** ROUTINE asm_inc_init *** … … 112 112 !! ** Action : 113 113 !!---------------------------------------------------------------------- 114 INTEGER, INTENT(in) :: Kmm ! time level index 114 INTEGER, INTENT(in) :: Kbb, Kmm, Krhs ! time level indices 115 ! 115 116 INTEGER :: ji, jj, jk, jt ! dummy loop indices 116 117 INTEGER :: imid, inum ! local integers … … 416 417 DO jj = 2, jpjm1 417 418 DO ji = fs_2, fs_jpim1 ! vector opt. 418 zhdiv(ji,jj) = ( e2u(ji ,jj) * e3u _n(ji ,jj,jk) * u_bkginc(ji ,jj,jk) &419 & - e2u(ji-1,jj) * e3u _n(ji-1,jj,jk) * u_bkginc(ji-1,jj,jk) &420 & + e1v(ji,jj ) * e3v _n(ji,jj ,jk) * v_bkginc(ji,jj ,jk) &421 & - e1v(ji,jj-1) * e3v _n(ji,jj-1,jk) * v_bkginc(ji,jj-1,jk) ) / e3t_n(ji,jj,jk)419 zhdiv(ji,jj) = ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * u_bkginc(ji ,jj,jk) & 420 & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * u_bkginc(ji-1,jj,jk) & 421 & + e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * v_bkginc(ji,jj ,jk) & 422 & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk) ) / e3t(ji,jj,jk,Kmm) 422 423 END DO 423 424 END DO … … 495 496 ! 496 497 IF( lk_asminc ) THEN !== data assimilation ==! 497 IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1 ) ! Output background fields498 IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1, Kmm ) ! Output background fields 498 499 IF( ln_asmdin ) THEN ! Direct initialization 499 IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1, K mm) ! Tracers500 IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 ) ! Dynamics501 IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 )! SSH500 IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1, Kbb, Kmm, ts , Krhs ) ! Tracers 501 IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1, Kbb, Kmm, uu, vv, Krhs ) ! Dynamics 502 IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1, Kbb, Kmm ) ! SSH 502 503 ENDIF 503 504 ENDIF … … 506 507 507 508 508 SUBROUTINE tra_asm_inc( kt, K mm)509 SUBROUTINE tra_asm_inc( kt, Kbb, Kmm, pts, Krhs ) 509 510 !!---------------------------------------------------------------------- 510 511 !! *** ROUTINE tra_asm_inc *** … … 516 517 !! ** Action : 517 518 !!---------------------------------------------------------------------- 518 INTEGER, INTENT(IN) :: kt ! Current time step 519 INTEGER, INTENT(IN) :: Kmm ! Current time level index 519 INTEGER , INTENT(in ) :: kt ! Current time step 520 INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! Time level indices 521 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 520 522 ! 521 523 INTEGER :: ji, jj, jk … … 528 530 ! used to prevent the applied increments taking the temperature below the local freezing point 529 531 DO jk = 1, jpkm1 530 CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), gdept_n(:,:,jk) )532 CALL eos_fzp( pts(:,:,jk,jp_sal,Kmm), fzptnz(:,:,jk), gdept(:,:,jk,Kmm) ) 531 533 END DO 532 534 ! … … 551 553 ! Do not apply negative increments if the temperature will fall below freezing 552 554 WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. & 553 & tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) )554 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt555 & pts(:,:,jk,jp_tem,Kmm) + pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 556 pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 555 557 END WHERE 556 558 ELSE 557 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt559 pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 558 560 ENDIF 559 561 IF (ln_salfix) THEN … … 561 563 ! minimum value salfixmin 562 564 WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. & 563 & tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin )564 tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt565 & pts(:,:,jk,jp_sal,Kmm) + pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 566 pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt 565 567 END WHERE 566 568 ELSE 567 tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt569 pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt 568 570 ENDIF 569 571 END DO … … 586 588 IF (ln_temnofreeze) THEN 587 589 ! Do not apply negative increments if the temperature will fall below freezing 588 WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) )589 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)590 WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_tem,Kmm) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 591 pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:) 590 592 END WHERE 591 593 ELSE 592 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)594 pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:) 593 595 ENDIF 594 596 IF (ln_salfix) THEN 595 597 ! Do not apply negative increments if the salinity will fall below a specified 596 598 ! minimum value salfixmin 597 WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin )598 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)599 WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_sal,Kmm) + s_bkginc(:,:,:) > salfixmin ) 600 pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:) 599 601 END WHERE 600 602 ELSE 601 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)603 pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:) 602 604 ENDIF 603 605 604 tsb(:,:,:,:) = tsn(:,:,:,:) ! Update before fields605 606 CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! Before potential and in situ densities606 pts(:,:,:,:,Kbb) = pts(:,:,:,:,Kmm) ! Update before fields 607 608 CALL eos( pts(:,:,:,:,Kbb), rhd, rhop, gdept_0(:,:,:) ) ! Before potential and in situ densities 607 609 !!gm fabien 608 ! CALL eos( tsb, rhd, rhop ) ! Before potential and in situ densities610 ! CALL eos( pts(:,:,:,:,Kbb), rhd, rhop ) ! Before potential and in situ densities 609 611 !!gm 610 612 611 613 IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav) & 612 & CALL zps_hde ( kt, Kmm, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient614 & CALL zps_hde ( kt, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, & ! Partial steps: before horizontal gradient 613 615 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 614 616 IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) & 615 & CALL zps_hde_isf( nit000, Kmm, jpts, tsb, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)617 & CALL zps_hde_isf( nit000, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 616 618 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the last ocean level 617 619 … … 629 631 630 632 631 SUBROUTINE dyn_asm_inc( kt )633 SUBROUTINE dyn_asm_inc( kt, Kbb, Kmm, puu, pvv, Krhs ) 632 634 !!---------------------------------------------------------------------- 633 635 !! *** ROUTINE dyn_asm_inc *** … … 639 641 !! ** Action : 640 642 !!---------------------------------------------------------------------- 641 INTEGER, INTENT(IN) :: kt ! Current time step 643 INTEGER , INTENT( in ) :: kt ! ocean time-step index 644 INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs ! ocean time level indices 645 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 642 646 ! 643 647 INTEGER :: jk … … 663 667 ! Update the dynamic tendencies 664 668 DO jk = 1, jpkm1 665 ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt666 va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt669 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + u_bkginc(:,:,jk) * zincwgt 670 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + v_bkginc(:,:,jk) * zincwgt 667 671 END DO 668 672 ! … … 682 686 ! 683 687 ! Initialize the now fields with the background + increment 684 un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:)685 vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:)686 ! 687 ub(:,:,:) = un(:,:,:) ! Update before fields688 vb(:,:,:) = vn(:,:,:)688 puu(:,:,:,Kmm) = u_bkg(:,:,:) + u_bkginc(:,:,:) 689 pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:) 690 ! 691 puu(:,:,:,Kbb) = puu(:,:,:,Kmm) ! Update before fields 692 pvv(:,:,:,Kbb) = pvv(:,:,:,Kmm) 689 693 ! 690 694 DEALLOCATE( u_bkg ) … … 699 703 700 704 701 SUBROUTINE ssh_asm_inc( kt )705 SUBROUTINE ssh_asm_inc( kt, Kbb, Kmm ) 702 706 !!---------------------------------------------------------------------- 703 707 !! *** ROUTINE ssh_asm_inc *** … … 709 713 !! ** Action : 710 714 !!---------------------------------------------------------------------- 711 INTEGER, INTENT(IN) :: kt ! Current time step 715 INTEGER, INTENT(IN) :: kt ! Current time step 716 INTEGER, INTENT(IN) :: Kbb, Kmm ! Current time step 712 717 ! 713 718 INTEGER :: it … … 756 761 neuler = 0 ! Force Euler forward step 757 762 ! 758 ssh n(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) ! Initialize the now fields the background + increment759 ! 760 ssh b(:,:) = sshn(:,:) ! Update before fields761 e3t _b(:,:,:) = e3t_n(:,:,:)762 !!gm why not e3u _b, e3v_b, gdept_b????763 ssh(:,:,Kmm) = ssh_bkg(:,:) + ssh_bkginc(:,:) ! Initialize the now fields the background + increment 764 ! 765 ssh(:,:,Kbb) = ssh(:,:,Kmm) ! Update before fields 766 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 767 !!gm why not e3u(:,:,:,Kbb), e3v(:,:,:,Kbb), gdept(:,:,:,Kbb) ???? 763 768 ! 764 769 DEALLOCATE( ssh_bkg ) … … 772 777 773 778 774 SUBROUTINE ssh_asm_div( kt, phdivn )779 SUBROUTINE ssh_asm_div( kt, Kbb, Kmm, phdivn ) 775 780 !!---------------------------------------------------------------------- 776 781 !! *** ROUTINE ssh_asm_div *** … … 786 791 !!---------------------------------------------------------------------- 787 792 INTEGER, INTENT(IN) :: kt ! ocean time-step index 793 INTEGER, INTENT(IN) :: Kbb, Kmm ! time level indices 788 794 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: phdivn ! horizontal divergence 789 795 !! … … 793 799 ! 794 800 #if defined key_asminc 795 CALL ssh_asm_inc( kt ) !== (calculate increments)801 CALL ssh_asm_inc( kt, Kbb, Kmm ) !== (calculate increments) 796 802 ! 797 803 IF( ln_linssh ) THEN 798 phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t _n(:,:,1) * tmask(:,:,1)804 phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t(:,:,1,Kmm) * tmask(:,:,1) 799 805 ELSE 800 806 ALLOCATE( ztim(jpi,jpj) )
Note: See TracChangeset
for help on using the changeset viewer.
