Changeset 15381
- Timestamp:
- 2021-10-15T12:03:54+02:00 (3 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE
- Files:
-
- 2 edited
-
DIA/diaregmean.F90 (modified) (1 diff)
-
ZDF/zdfmxl.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/DIA/diaregmean.F90
r15378 r15381 648 648 649 649 650 !JT MLDIF( ln_diaregmean_karamld ) THEN651 !JT MLDtmp_field_mat(:,:,16) = tmp_field_mat(:,:,16) + (hmld_kara(:,:)*tmask(:,:,1)) !mldkara652 !JT MLDENDIF650 IF( ln_diaregmean_karamld ) THEN 651 tmp_field_mat(:,:,16) = tmp_field_mat(:,:,16) + (hmld_kara(:,:)*tmask(:,:,1)) !mldkara 652 ENDIF 653 653 654 654 name_dat_mat(16) = 'mldkara' -
NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/ZDF/zdfmxl.F90
r14078 r15381 22 22 USE iom ! I/O library 23 23 USE lib_mpp ! MPP library 24 !JT 25 USE lbclnk ! (or mpp link) 26 !JT 24 27 25 28 IMPLICIT NONE … … 36 39 LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ll_found ! Is T_b to be found by interpolation ? 37 40 LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ll_belowml ! Flag points below mixed layer when ll_found=F 41 42 !JT 43 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmld_kara !: mixed layer depth of Kara et al [m] 44 45 46 47 ! Namelist variables for namzdf_karaml 48 LOGICAL :: ln_kara ! Logical switch for calculating kara mixed 49 ! layer 50 LOGICAL :: ln_kara_write25h ! Logical switch for 25 hour outputs 51 INTEGER :: jpmld_type ! mixed layer type 52 REAL(wp) :: ppz_ref ! depth of initial T_ref 53 REAL(wp) :: ppdT_crit ! Critical temp diff 54 REAL(wp) :: ppiso_frac ! Fraction of ppdT_crit used 55 56 !Used for 25h mean 57 LOGICAL, PRIVATE :: kara_25h_init = .TRUE. !Logical used to initalise 25h 58 !outputs. Necissary, because we need to 59 !initalise the kara_25h on the zeroth 60 !timestep (i.e in the nemogcm_init call) 61 REAL, PRIVATE, ALLOCATABLE, DIMENSION(:,:) :: hmld_kara_25h 62 63 64 65 !JT 38 66 39 67 REAL(wp), PUBLIC :: rho_c = 0.01_wp !: density criterion for mixed layer depth … … 127 155 END DO 128 156 ! depth of the mixing and mixed layers 157 !JT 158 CALL zdf_mxl_kara( kt ) ! kara MLD 159 !JT 160 129 161 DO jj = 1, jpj 130 162 DO ji = 1, jpi … … 489 521 END SUBROUTINE zdf_mxl_zint 490 522 523 524 525 SUBROUTINE zdf_mxl_kara( kt ) 526 !!---------------------------------------------------------------------------------- 527 !! *** ROUTINE zdf_mxl_kara *** 528 ! 529 ! Calculate mixed layer depth according to the definition of 530 ! Kara et al, 2000, JGR, 105, 16803. 531 ! 532 ! If mld_type=1 the mixed layer depth is calculated as the depth at which the 533 ! density has increased by an amount equivalent to a temperature difference of 534 ! 0.8C at the surface. 535 ! 536 ! For other values of mld_type the mixed layer is calculated as the depth at 537 ! which the temperature differs by 0.8C from the surface temperature. 538 ! 539 ! Original version: David Acreman 540 ! 541 !!----------------------------------------------------------------------------------- 542 543 INTEGER, INTENT( in ) :: kt ! ocean time-step index 544 545 NAMELIST/namzdf_karaml/ ln_kara,jpmld_type, ppz_ref, ppdT_crit, & 546 & ppiso_frac, ln_kara_write25h 547 548 ! Local variables 549 REAL, DIMENSION(jpi,jpk) :: ppzdep ! depth for use in calculating d(rho) 550 REAL(wp), DIMENSION(jpi,jpj,jpts) :: ztsn_2d !Local version of tsn 551 LOGICAL :: ll_found(jpi,jpj) ! Is T_b to be found by interpolation ? 552 LOGICAL :: ll_belowml(jpi,jpj,jpk) ! Flag points below mixed layer when ll_found=F 553 INTEGER :: ji, jj, jk ! loop counter 554 INTEGER :: ik_ref(jpi,jpj) ! index of reference level 555 INTEGER :: ik_iso(jpi,jpj) ! index of last uniform temp level 556 REAL :: zT(jpi,jpj,jpk) ! Temperature or denisty 557 REAL :: zT_ref(jpi,jpj) ! reference temperature 558 REAL :: zT_b ! base temperature 559 REAL :: zdTdz(jpi,jpj,jpk-2) ! gradient of zT 560 REAL :: zmoddT(jpi,jpj,jpk-2) ! Absolute temperature difference 561 REAL :: zdz ! depth difference 562 REAL :: zdT ! temperature difference 563 REAL :: zdelta_T(jpi,jpj) ! difference critereon 564 REAL :: zRHO1(jpi,jpj), zRHO2(jpi,jpj) ! Densities 565 INTEGER, SAVE :: idt, i_steps 566 INTEGER, SAVE :: i_cnt_25h ! Counter for 25 hour means 567 INTEGER :: ios ! Local integer output status for namelist read 568 569 570 571 !!------------------------------------------------------------------------------------- 572 573 IF( kt == nit000 ) THEN 574 ! Default values 575 ln_kara = .FALSE. 576 ln_kara_write25h = .FALSE. 577 jpmld_type = 1 578 ppz_ref = 10.0 579 ppdT_crit = 0.2 580 ppiso_frac = 0.1 581 ! Read namelist 582 REWIND ( numnam_ref ) 583 READ ( numnam_ref, namzdf_karaml, IOSTAT=ios, ERR= 901 ) 584 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_karaml in reference namelist' ) 585 586 REWIND( numnam_cfg ) ! Namelist nam_diadiaop in configuration namelist 3D hourly diagnostics 587 READ ( numnam_cfg, namzdf_karaml, IOSTAT = ios, ERR = 902 ) 588 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_karaml in configuration namelist' ) 589 IF(lwm) WRITE ( numond, namzdf_karaml ) 590 591 592 WRITE(numout,*) '===== Kara mixed layer =====' 593 WRITE(numout,*) 'ln_kara = ', ln_kara 594 WRITE(numout,*) 'jpmld_type = ', jpmld_type 595 WRITE(numout,*) 'ppz_ref = ', ppz_ref 596 WRITE(numout,*) 'ppdT_crit = ', ppdT_crit 597 WRITE(numout,*) 'ppiso_frac = ', ppiso_frac 598 WRITE(numout,*) 'ln_kara_write25h = ', ln_kara_write25h 599 WRITE(numout,*) '============================' 600 601 IF ( .NOT.ln_kara ) THEN 602 WRITE(numout,*) "ln_kara not set; Kara mixed layer not calculated" 603 ELSE 604 IF (.NOT. ALLOCATED(hmld_kara) ) ALLOCATE( hmld_kara(jpi,jpj) ) 605 IF ( ln_kara_write25h .AND. kara_25h_init ) THEN 606 i_cnt_25h = 0 607 IF (.NOT. ALLOCATED(hmld_kara_25h) ) & 608 & ALLOCATE( hmld_kara_25h(jpi,jpj) ) 609 hmld_kara_25h = 0._wp 610 !IF( nacc == 1 ) THEN 611 ! idt = INT(rdtmin) 612 !ELSE 613 ! idt = INT(rdt) 614 !ENDIF 615 616 idt = INT(rdt) 617 IF( MOD( 3600,idt) == 0 ) THEN 618 i_steps = 3600 / idt 619 ELSE 620 CALL ctl_stop('STOP', & 621 & 'zdf_mxl_kara: timestep must give MOD(3600,rdt)'// & 622 & ' = 0 otherwise no hourly values are possible') 623 ENDIF 624 ENDIF 625 ENDIF 626 ENDIF 627 628 IF ( ln_kara ) THEN 629 630 !set critical ln_kara 631 ztsn_2d = tsn(:,:,1,:) 632 ztsn_2d(:,:,jp_tem) = ztsn_2d(:,:,jp_tem) + ppdT_crit 633 634 ! Set the mixed layer depth criterion at each grid point 635 ppzdep = 0._wp 636 IF( jpmld_type == 1 ) THEN 637 CALL eos ( tsn(:,:,1,:), & 638 & ppzdep(:,:), zRHO1(:,:) ) 639 CALL eos ( ztsn_2d(:,:,:), & 640 & ppzdep(:,:), zRHO2(:,:) ) 641 zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0 642 ! RHO from eos (2d version) doesn't calculate north or east halo: 643 CALL lbc_lnk( 'zdf_mxl_kara',zdelta_T, 'T', 1. ) 644 zT(:,:,:) = rhop(:,:,:) 645 ELSE 646 zdelta_T(:,:) = ppdT_crit 647 zT(:,:,:) = tsn(:,:,:,jp_tem) 648 ENDIF 649 650 ! Calculate the gradient of zT and absolute difference for use later 651 DO jk = 1 ,jpk-2 652 zdTdz(:,:,jk) = ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1) 653 zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) ) 654 END DO 655 656 ! Find density/temperature at the reference level (Kara et al use 10m). 657 ! ik_ref is the index of the box centre immediately above or at the reference level 658 ! Find ppz_ref in the array of model level depths and find the ref 659 ! density/temperature by linear interpolation. 660 ik_ref = -1 661 DO jk = jpkm1, 2, -1 662 WHERE( gdept_n(:,:,jk) > ppz_ref ) 663 ik_ref(:,:) = jk - 1 664 zT_ref(:,:) = zT(:,:,jk-1) + & 665 & zdTdz(:,:,jk-1) * ( ppz_ref - gdept_n(:,:,jk-1) ) 666 ENDWHERE 667 END DO 668 IF ( ANY( ik_ref < 0 ) .OR. ANY( ik_ref > jpkm1 ) ) THEN 669 CALL ctl_stop( "STOP", & 670 & "zdf_mxl_kara: unable to find reference level for kara ML" ) 671 ELSE 672 ! If the first grid box centre is below the reference level then use the 673 ! top model level to get zT_ref 674 WHERE( gdept_n(:,:,1) > ppz_ref ) 675 zT_ref = zT(:,:,1) 676 ik_ref = 1 677 ENDWHERE 678 679 ! Search for a uniform density/temperature region where adjacent levels 680 ! differ by less than ppiso_frac * deltaT. 681 ! ik_iso is the index of the last level in the uniform layer 682 ! ll_found indicates whether the mixed layer depth can be found by interpolation 683 ik_iso(:,:) = ik_ref(:,:) 684 ll_found(:,:) = .false. 685 DO jj = 1, nlcj 686 DO ji = 1, nlci 687 !CDIR NOVECTOR 688 DO jk = ik_ref(ji,jj), mbkt(ji,jj)-1 !mbathy(ji,jj)-1 689 IF( zmoddT(ji,jj,jk) > ( ppiso_frac * zdelta_T(ji,jj) ) ) THEN 690 ik_iso(ji,jj) = jk 691 ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 692 EXIT 693 ENDIF 694 END DO 695 END DO 696 END DO 697 698 ! Use linear interpolation to find depth of mixed layer base where possible 699 hmld_kara(:,:) = ppz_ref 700 DO jj = 1, jpj 701 DO ji = 1, jpi 702 IF( ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0 ) THEN 703 zdz = abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 704 hmld_kara(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz 705 ENDIF 706 END DO 707 END DO 708 709 ! If ll_found = .false. then calculate MLD using difference of zdelta_T 710 ! from the reference density/temperature 711 712 ! Prevent this section from working on land points 713 WHERE( tmask(:,:,1) /= 1.0 ) 714 ll_found = .true. 715 ENDWHERE 716 717 DO jk = 1, jpk 718 ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= & 719 & zdelta_T(:,:) 720 END DO 721 722 ! Set default value where interpolation cannot be used (ll_found=false) 723 DO jj = 1, jpj 724 DO ji = 1, jpi 725 IF( .NOT. ll_found(ji,jj) ) & 726 & hmld_kara(ji,jj) = gdept_n(ji,jj,mbkt(ji,jj))! mbathy(ji,jj)) 727 END DO 728 END DO 729 730 DO jj = 1, jpj 731 DO ji = 1, jpi 732 !CDIR NOVECTOR 733 DO jk = ik_ref(ji,jj)+1, mbkt(ji,jj) !mbathy(ji,jj) 734 IF( ll_found(ji,jj) ) EXIT 735 IF( ll_belowml(ji,jj,jk) ) THEN 736 zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * & 737 & SIGN(1.0, zdTdz(ji,jj,jk-1) ) 738 zdT = zT_b - zT(ji,jj,jk-1) 739 zdz = zdT / zdTdz(ji,jj,jk-1) 740 hmld_kara(ji,jj) = gdept_n(ji,jj,jk-1) + zdz 741 EXIT 742 ENDIF 743 END DO 744 END DO 745 END DO 746 747 hmld_kara(:,:) = hmld_kara(:,:) * tmask(:,:,1) 748 749 IF( ln_kara_write25h ) THEN 750 !Double IF required as i_steps not defined if ln_kara_write25h = 751 ! FALSE 752 IF ( ( MOD( kt, i_steps ) == 0 ) .OR. kara_25h_init ) THEN 753 hmld_kara_25h = hmld_kara_25h + hmld_kara 754 i_cnt_25h = i_cnt_25h + 1 755 IF ( kara_25h_init ) kara_25h_init = .FALSE. 756 ENDIF 757 ENDIF 758 759 !#if defined key_iomput 760 IF( kt /= nit000 ) THEN 761 CALL iom_put( "mldkara" , hmld_kara ) 762 IF( ( MOD( i_cnt_25h, 25) == 0 ) .AND. ln_kara_write25h ) & 763 CALL iom_put( "kara25h" , ( hmld_kara_25h / 25._wp ) ) 764 ENDIF 765 !#endif 766 767 ENDIF 768 ENDIF 769 770 END SUBROUTINE zdf_mxl_kara 771 491 772 !!====================================================================== 492 773 END MODULE zdfmxl
Note: See TracChangeset
for help on using the changeset viewer.
