- Timestamp:
- 2020-06-24T14:38:26+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/ASM/asminc.F90
r12489 r13151 9 9 !! ! 2007-04 (A. Weaver) Merge with OPAVAR/NEMOVAR 10 10 !! NEMO 3.3 ! 2010-05 (D. Lea) Update to work with NEMO v3.2 11 !! - ! 2010-05 (D. Lea) add calc_month_len routine based on day_init 11 !! - ! 2010-05 (D. Lea) add calc_month_len routine based on day_init 12 12 !! 3.4 ! 2012-10 (A. Weaver and K. Mogensen) Fix for direct initialization 13 13 !! ! 2014-09 (D. Lea) Local calc_date removed use routine from OBS … … 31 31 USE zpshde ! Partial step : Horizontal Derivative 32 32 USE asmpar ! Parameters for the assmilation interface 33 USE asmbkg ! 33 USE asmbkg ! 34 34 USE c1d ! 1D initialization 35 35 USE sbc_oce ! Surface boundary condition variables. … … 45 45 IMPLICIT NONE 46 46 PRIVATE 47 47 48 48 PUBLIC asm_inc_init !: Initialize the increment arrays and IAU weights 49 49 PUBLIC tra_asm_inc !: Apply the tracer (T and S) increments … … 72 72 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkg , v_bkg !: Background u- & v- velocity components 73 73 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkginc, s_bkginc !: Increment to the background T & S 74 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u- & v-components 74 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u- & v-components 75 75 REAL(wp), PUBLIC, DIMENSION(:) , ALLOCATABLE :: wgtiau !: IAU weights for each time step 76 76 #if defined key_asminc … … 80 80 INTEGER , PUBLIC :: nitbkg !: Time step of the background state used in the Jb term 81 81 INTEGER , PUBLIC :: nitdin !: Time step of the background state for direct initialization 82 INTEGER , PUBLIC :: nitiaustr !: Time step of the start of the IAU interval 82 INTEGER , PUBLIC :: nitiaustr !: Time step of the start of the IAU interval 83 83 INTEGER , PUBLIC :: nitiaufin !: Time step of the end of the IAU interval 84 ! 84 ! 85 85 INTEGER , PUBLIC :: niaufn !: Type of IAU weighing function: = 0 Constant weighting 86 ! !: = 1 Linear hat-like, centred in middle of IAU interval 86 ! !: = 1 Linear hat-like, centred in middle of IAU interval 87 87 REAL(wp), PUBLIC :: salfixmin !: Ensure that the salinity is larger than this value if (ln_salfix) 88 88 … … 95 95 !! * Substitutions 96 96 # include "do_loop_substitute.h90" 97 # include "domzgr_substitute.h90" 97 98 !!---------------------------------------------------------------------- 98 99 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 105 106 !!---------------------------------------------------------------------- 106 107 !! *** ROUTINE asm_inc_init *** 107 !! 108 !! 108 109 !! ** Purpose : Initialize the assimilation increment and IAU weights. 109 110 !! 110 111 !! ** Method : Initialize the assimilation increment and IAU weights. 111 112 !! 112 !! ** Action : 113 !! ** Action : 113 114 !!---------------------------------------------------------------------- 114 115 INTEGER, INTENT(in) :: Kbb, Kmm, Krhs ! time level indices … … 262 263 ! 263 264 ! !--------------------------------------------------------- 264 IF( niaufn == 0 ) THEN ! Constant IAU forcing 265 IF( niaufn == 0 ) THEN ! Constant IAU forcing 265 266 ! !--------------------------------------------------------- 266 267 DO jt = 1, iiauper … … 268 269 END DO 269 270 ! !--------------------------------------------------------- 270 ELSEIF ( niaufn == 1 ) THEN ! Linear hat-like, centred in middle of IAU interval 271 ELSEIF ( niaufn == 1 ) THEN ! Linear hat-like, centred in middle of IAU interval 271 272 ! !--------------------------------------------------------- 272 273 ! Compute the normalization factor 273 274 znorm = 0._wp 274 275 IF( MOD( iiauper, 2 ) == 0 ) THEN ! Even number of time steps in IAU interval 275 imid = iiauper / 2 276 imid = iiauper / 2 276 277 DO jt = 1, imid 277 278 znorm = znorm + REAL( jt ) … … 279 280 znorm = 2.0 * znorm 280 281 ELSE ! Odd number of time steps in IAU interval 281 imid = ( iiauper + 1 ) / 2 282 imid = ( iiauper + 1 ) / 2 282 283 DO jt = 1, imid - 1 283 284 znorm = znorm + REAL( jt ) … … 306 307 DO jt = 1, icycper 307 308 ztotwgt = ztotwgt + wgtiau(jt) 308 WRITE(numout,*) ' ', jt, ' ', wgtiau(jt) 309 END DO 309 WRITE(numout,*) ' ', jt, ' ', wgtiau(jt) 310 END DO 310 311 WRITE(numout,*) ' ===================================' 311 312 WRITE(numout,*) ' Time-integrated weight = ', ztotwgt 312 313 WRITE(numout,*) ' ===================================' 313 314 ENDIF 314 315 315 316 ENDIF 316 317 … … 337 338 CALL iom_open( c_asminc, inum ) 338 339 ! 339 CALL iom_get( inum, 'time' , zdate_inc ) 340 CALL iom_get( inum, 'time' , zdate_inc ) 340 341 CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb ) 341 342 CALL iom_get( inum, 'z_inc_datef', z_inc_datef ) … … 344 345 ! 345 346 IF(lwp) THEN 346 WRITE(numout,*) 347 WRITE(numout,*) 347 348 WRITE(numout,*) 'asm_inc_init : Assimilation increments valid between dates ', z_inc_dateb,' and ', z_inc_datef 348 349 WRITE(numout,*) '~~~~~~~~~~~~' … … 358 359 & ' not agree with Direct Initialization time' ) 359 360 360 IF ( ln_trainc ) THEN 361 IF ( ln_trainc ) THEN 361 362 CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 ) 362 363 CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 ) … … 370 371 ENDIF 371 372 372 IF ( ln_dyninc ) THEN 373 CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 ) 374 CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 ) 373 IF ( ln_dyninc ) THEN 374 CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 ) 375 CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 ) 375 376 ! Apply the masks 376 377 u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:) … … 381 382 WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0 382 383 ENDIF 383 384 384 385 IF ( ln_sshinc ) THEN 385 386 CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 ) … … 407 408 IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN ! Apply divergence damping filter 408 409 ! !-------------------------------------- 409 ALLOCATE( zhdiv(jpi,jpj) ) 410 ALLOCATE( zhdiv(jpi,jpj) ) 410 411 ! 411 412 DO jt = 1, nn_divdmp … … 417 418 & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * u_bkginc(ji-1,jj,jk) & 418 419 & + e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * v_bkginc(ji,jj ,jk) & 419 & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk) ) / e3t(ji,jj,jk,Kmm) 420 & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk) ) & 421 & / e3t(ji,jj,jk,Kmm) 420 422 END_2D 421 423 CALL lbc_lnk( 'asminc', zhdiv, 'T', 1. ) ! lateral boundary cond. (no sign change) … … 425 427 & + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 426 428 v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) & 427 & + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 429 & + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 428 430 END_2D 429 431 END DO … … 431 433 END DO 432 434 ! 433 DEALLOCATE( zhdiv ) 435 DEALLOCATE( zhdiv ) 434 436 ! 435 437 ENDIF … … 452 454 CALL iom_open( c_asmdin, inum ) 453 455 ! 454 CALL iom_get( inum, 'rdastp', zdate_bkg ) 456 CALL iom_get( inum, 'rdastp', zdate_bkg ) 455 457 ! 456 458 IF(lwp) THEN 457 WRITE(numout,*) 459 WRITE(numout,*) 458 460 WRITE(numout,*) ' ==>>> Assimilation background state valid at : ', zdate_bkg 459 461 WRITE(numout,*) … … 464 466 & ' not agree with Direct Initialization time' ) 465 467 ! 466 IF ( ln_trainc ) THEN 468 IF ( ln_trainc ) THEN 467 469 CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg ) 468 470 CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg ) … … 471 473 ENDIF 472 474 ! 473 IF ( ln_dyninc ) THEN 475 IF ( ln_dyninc ) THEN 474 476 CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg ) 475 477 CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg ) … … 499 501 ! 500 502 END SUBROUTINE asm_inc_init 501 502 503 504 503 505 SUBROUTINE tra_asm_inc( kt, Kbb, Kmm, pts, Krhs ) 504 506 !!---------------------------------------------------------------------- 505 507 !! *** ROUTINE tra_asm_inc *** 506 !! 508 !! 507 509 !! ** Purpose : Apply the tracer (T and S) assimilation increments 508 510 !! 509 511 !! ** Method : Direct initialization or Incremental Analysis Updating 510 512 !! 511 !! ** Action : 513 !! ** Action : 512 514 !!---------------------------------------------------------------------- 513 515 INTEGER , INTENT(in ) :: kt ! Current time step … … 521 523 !!---------------------------------------------------------------------- 522 524 ! 523 ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 524 ! used to prevent the applied increments taking the temperature below the local freezing point 525 ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 526 ! used to prevent the applied increments taking the temperature below the local freezing point 525 527 DO jk = 1, jpkm1 526 528 CALL eos_fzp( pts(:,:,jk,jp_sal,Kmm), fzptnz(:,:,jk), gdept(:,:,jk,Kmm) ) … … 537 539 ! 538 540 IF(lwp) THEN 539 WRITE(numout,*) 541 WRITE(numout,*) 540 542 WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 541 543 WRITE(numout,*) '~~~~~~~~~~~~' … … 547 549 ! Do not apply negative increments if the temperature will fall below freezing 548 550 WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. & 549 & pts(:,:,jk,jp_tem,Kmm) + pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 550 pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 551 & pts(:,:,jk,jp_tem,Kmm) + pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 552 pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 551 553 END WHERE 552 554 ELSE 553 pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 555 pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 554 556 ENDIF 555 557 IF (ln_salfix) THEN … … 557 559 ! minimum value salfixmin 558 560 WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. & 559 & pts(:,:,jk,jp_sal,Kmm) + pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 561 & pts(:,:,jk,jp_sal,Kmm) + pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 560 562 pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt 561 563 END WHERE … … 574 576 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 575 577 ! !-------------------------------------- 576 ! 578 ! 577 579 IF ( kt == nitdin_r ) THEN 578 580 ! … … 582 584 IF (ln_temnofreeze) THEN 583 585 ! Do not apply negative increments if the temperature will fall below freezing 584 WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_tem,Kmm) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 585 pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:) 586 WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_tem,Kmm) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 587 pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:) 586 588 END WHERE 587 589 ELSE 588 pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:) 590 pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:) 589 591 ENDIF 590 592 IF (ln_salfix) THEN 591 593 ! Do not apply negative increments if the salinity will fall below a specified 592 594 ! minimum value salfixmin 593 WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_sal,Kmm) + s_bkginc(:,:,:) > salfixmin ) 594 pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:) 595 WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_sal,Kmm) + s_bkginc(:,:,:) > salfixmin ) 596 pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:) 595 597 END WHERE 596 598 ELSE 597 pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:) 599 pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:) 598 600 ENDIF 599 601 … … 617 619 DEALLOCATE( s_bkg ) 618 620 ENDIF 619 ! 621 ! 620 622 ENDIF 621 623 ! Perhaps the following call should be in step … … 628 630 !!---------------------------------------------------------------------- 629 631 !! *** ROUTINE dyn_asm_inc *** 630 !! 632 !! 631 633 !! ** Purpose : Apply the dynamics (u and v) assimilation increments. 632 634 !! 633 635 !! ** Method : Direct initialization or Incremental Analysis Updating. 634 636 !! 635 !! ** Action : 637 !! ** Action : 636 638 !!---------------------------------------------------------------------- 637 639 INTEGER , INTENT( in ) :: kt ! ocean time-step index … … 654 656 ! 655 657 IF(lwp) THEN 656 WRITE(numout,*) 658 WRITE(numout,*) 657 659 WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 658 660 WRITE(numout,*) '~~~~~~~~~~~~' … … 674 676 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 675 677 ! !----------------------------------------- 676 ! 678 ! 677 679 IF ( kt == nitdin_r ) THEN 678 680 ! … … 681 683 ! Initialize the now fields with the background + increment 682 684 puu(:,:,:,Kmm) = u_bkg(:,:,:) + u_bkginc(:,:,:) 683 pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:) 685 pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:) 684 686 ! 685 687 puu(:,:,:,Kbb) = puu(:,:,:,Kmm) ! Update before fields … … 700 702 !!---------------------------------------------------------------------- 701 703 !! *** ROUTINE ssh_asm_inc *** 702 !! 704 !! 703 705 !! ** Purpose : Apply the sea surface height assimilation increment. 704 706 !! 705 707 !! ** Method : Direct initialization or Incremental Analysis Updating. 706 708 !! 707 !! ** Action : 709 !! ** Action : 708 710 !!---------------------------------------------------------------------- 709 711 INTEGER, INTENT(IN) :: kt ! Current time step … … 725 727 ! 726 728 IF(lwp) THEN 727 WRITE(numout,*) 729 WRITE(numout,*) 728 730 WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', & 729 731 & kt,' with IAU weight = ', wgtiau(it) … … 758 760 ! 759 761 ssh(:,:,Kbb) = ssh(:,:,Kmm) ! Update before fields 762 #if ! defined key_qco 760 763 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 764 #endif 761 765 !!gm why not e3u(:,:,:,Kbb), e3v(:,:,:,Kbb), gdept(:,:,:,Kbb) ???? 762 766 ! … … 775 779 !! *** ROUTINE ssh_asm_div *** 776 780 !! 777 !! ** Purpose : ssh increment with z* is incorporated via a correction of the local divergence 781 !! ** Purpose : ssh increment with z* is incorporated via a correction of the local divergence 778 782 !! across all the water column 779 783 !! … … 791 795 REAL(wp), DIMENSION(:,:) , POINTER :: ztim ! local array 792 796 !!---------------------------------------------------------------------- 793 ! 797 ! 794 798 #if defined key_asminc 795 799 CALL ssh_asm_inc( kt, Kbb, Kmm ) !== (calculate increments) 796 800 ! 797 IF( ln_linssh ) THEN 801 IF( ln_linssh ) THEN 798 802 phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t(:,:,1,Kmm) * tmask(:,:,1) 799 ELSE 803 ELSE 800 804 ALLOCATE( ztim(jpi,jpj) ) 801 805 ztim(:,:) = ssh_iau(:,:) / ( ht(:,:) + 1.0 - ssmask(:,:) ) 802 DO jk = 1, jpkm1 803 phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) 806 DO jk = 1, jpkm1 807 phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) 804 808 END DO 805 809 ! … … 814 818 !!---------------------------------------------------------------------- 815 819 !! *** ROUTINE seaice_asm_inc *** 816 !! 820 !! 817 821 !! ** Purpose : Apply the sea ice assimilation increment. 818 822 !! 819 823 !! ** Method : Direct initialization or Incremental Analysis Updating. 820 824 !! 821 !! ** Action : 825 !! ** Action : 822 826 !! 823 827 !!---------------------------------------------------------------------- … … 840 844 ! 841 845 it = kt - nit000 + 1 842 zincwgt = wgtiau(it) ! IAU weight for the current time step 846 zincwgt = wgtiau(it) ! IAU weight for the current time step 843 847 ! note this is not a tendency so should not be divided by rn_Dt (as with the tracer and other increments) 844 848 ! 845 849 IF(lwp) THEN 846 WRITE(numout,*) 850 WRITE(numout,*) 847 851 WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 848 852 WRITE(numout,*) '~~~~~~~~~~~~' … … 862 866 ! 863 867 ! Nudge sea ice depth to bring it up to a required minimum depth 864 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 865 zhicifinc(:,:) = (zhicifmin - hm_i(:,:)) * zincwgt 868 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 869 zhicifinc(:,:) = (zhicifmin - hm_i(:,:)) * zincwgt 866 870 ELSEWHERE 867 871 zhicifinc(:,:) = 0.0_wp … … 896 900 IF ( kt == nitdin_r ) THEN 897 901 ! 902 <<<<<<< .working 898 903 l_1st_euler = 0 ! Force Euler forward step 904 ======= 905 l_1st_euler = .TRUE. ! Force Euler forward step 906 >>>>>>> .merge-right.r13092 899 907 ! 900 908 ! Sea-ice : SI3 case … … 903 911 zofrld (:,:) = 1._wp - at_i(:,:) 904 912 zohicif(:,:) = hm_i(:,:) 905 ! 913 ! 906 914 ! Initialize the now fields the background + increment 907 915 at_i(:,:) = 1. - MIN( MAX( 1.-at_i(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp) 908 at_i_b(:,:) = at_i(:,:) 916 at_i_b(:,:) = at_i(:,:) 909 917 fr_i(:,:) = at_i(:,:) ! adjust ice fraction 910 918 ! … … 912 920 ! 913 921 ! Nudge sea ice depth to bring it up to a required minimum depth 914 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 922 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 915 923 zhicifinc(:,:) = zhicifmin - hm_i(:,:) 916 924 ELSEWHERE … … 942 950 !#if defined defined key_si3 || defined key_cice 943 951 ! 944 ! IF (ln_seaicebal ) THEN 952 ! IF (ln_seaicebal ) THEN 945 953 ! !! balancing salinity increments 946 954 ! !! simple case from limflx.F90 (doesn't include a mass flux) … … 954 962 ! 955 963 ! DO jj = 1, jpj 956 ! DO ji = 1, jpi 964 ! DO ji = 1, jpi 957 965 ! ! calculate change in ice and snow mass per unit area 958 966 ! ! positive values imply adding salt to the ocean (results from ice formation) … … 965 973 ! 966 974 ! ! prevent small mld 967 ! ! less than 10m can cause salinity instability 975 ! ! less than 10m can cause salinity instability 968 976 ! IF (mld < 10) mld=10 969 977 ! 970 ! ! set to bottom of a level 978 ! ! set to bottom of a level 971 979 ! DO jk = jpk-1, 2, -1 972 ! IF ((mld > gdepw(ji,jj,jk )) .and. (mld < gdepw(ji,jj,jk+1))) THEN973 ! mld=gdepw(ji,jj,jk+1 )980 ! IF ((mld > gdepw(ji,jj,jk,Kmm)) .and. (mld < gdepw(ji,jj,jk+1,Kmm))) THEN 981 ! mld=gdepw(ji,jj,jk+1,Kmm) 974 982 ! jkmax=jk 975 983 ! ENDIF … … 977 985 ! 978 986 ! ! avoid applying salinity balancing in shallow water or on land 979 ! ! 987 ! ! 980 988 ! 981 989 ! ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m) … … 988 996 ! 989 997 ! ! put increments in for levels in the mixed layer 990 ! ! but prevent salinity below a threshold value 991 ! 992 ! DO jk = 1, jkmax 993 ! 994 ! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 998 ! ! but prevent salinity below a threshold value 999 ! 1000 ! DO jk = 1, jkmax 1001 ! 1002 ! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 995 1003 ! sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn 996 1004 ! sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn … … 1003 1011 ! ! 1004 1012 ! !! Adjust fsalt. A +ve fsalt means adding salt to ocean 1005 ! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt 1006 ! !! 1007 ! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d) 1013 ! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt 1014 ! !! 1015 ! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d) 1008 1016 ! !! ! E-P (kg m-2 s-2) 1009 1017 ! ! emp(ji,jj) = emp(ji,jj) + zpmess ! E-P (kg m-2 s-2) … … 1018 1026 ! 1019 1027 END SUBROUTINE seaice_asm_inc 1020 1028 1021 1029 !!====================================================================== 1022 1030 END MODULE asminc
Note: See TracChangeset
for help on using the changeset viewer.
