240
241
242
243
244 integer, intent(in) :: ng, tile
245 integer, intent(in) :: LBi, UBi, LBj, UBj, LBij, UBij
246 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
247 integer, intent(in) :: Lini
248
249 real(r8), intent(in) :: state(Ninner)
250
251# ifdef ASSUMED_SHAPE
252# ifdef MASKING
253 real(r8), intent(in) :: rmask(LBi:,LBj:)
254 real(r8), intent(in) :: umask(LBi:,LBj:)
255 real(r8), intent(in) :: vmask(LBi:,LBj:)
256# endif
257# ifdef ADJUST_BOUNDARY
258# ifdef SOLVE3D
259 real(r8), intent(inout) :: tl_t_obc(LBij:,:,:,:,:,:)
260 real(r8), intent(inout) :: tl_u_obc(LBij:,:,:,:,:)
261 real(r8), intent(inout) :: tl_v_obc(LBij:,:,:,:,:)
262# ifdef LCZ_FINAL
263 real(r8), intent(inout) :: ad_t_obc(LBij:,:,:,:,:,:)
264 real(r8), intent(inout) :: ad_u_obc(LBij:,:,:,:,:)
265 real(r8), intent(inout) :: ad_v_obc(LBij:,:,:,:,:)
266# endif
267# endif
268 real(r8), intent(inout) :: tl_ubar_obc(LBij:,:,:,:)
269 real(r8), intent(inout) :: tl_vbar_obc(LBij:,:,:,:)
270 real(r8), intent(inout) :: tl_zeta_obc(LBij:,:,:,:)
271# ifdef LCZ_FINAL
272 real(r8), intent(inout) :: ad_ubar_obc(LBij:,:,:,:)
273 real(r8), intent(inout) :: ad_vbar_obc(LBij:,:,:,:)
274 real(r8), intent(inout) :: ad_zeta_obc(LBij:,:,:,:)
275# endif
276# endif
277# ifdef ADJUST_WSTRESS
278 real(r8), intent(inout) :: tl_ustr(LBi:,LBj:,:,:)
279 real(r8), intent(inout) :: tl_vstr(LBi:,LBj:,:,:)
280# ifdef LCZ_FINAL
281 real(r8), intent(inout) :: ad_ustr(LBi:,LBj:,:,:)
282 real(r8), intent(inout) :: ad_vstr(LBi:,LBj:,:,:)
283# endif
284# endif
285# if defined ADJUST_STFLUX && defined SOLVE3D
286 real(r8), intent(inout) :: tl_tflux(LBi:,LBj:,:,:,:)
287# ifdef LCZ_FINAL
288 real(r8), intent(inout) :: ad_tflux(LBi:,LBj:,:,:,:)
289# endif
290# endif
291# ifdef SOLVE3D
292 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
293 real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:)
294 real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:)
295# if defined WEAK_CONSTRAINT && defined TIME_CONV
296 real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
297 real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
298# endif
299# ifdef LCZ_FINAL
300 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
301 real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:)
302 real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:)
303# if defined WEAK_CONSTRAINT && defined TIME_CONV
304 real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
305 real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
306# endif
307# endif
308# ifdef HESSIAN_FSV
309 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
310 real(r8), intent(inout) :: f_t(LBi:,LBj:,:,:)
311 real(r8), intent(inout) :: f_u(LBi:,LBj:,:)
312 real(r8), intent(inout) :: f_v(LBi:,LBj:,:)
313 real(r8), intent(inout) :: f_ubar(LBi:,LBj:)
314 real(r8), intent(inout) :: f_vbar(LBi:,LBj:)
315# endif
316# else
317 real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
318 real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
319# ifdef HESSIAN_FSV
320 real(r8), intent(inout) :: f_ubar(LBi:,LBj:)
321 real(r8), intent(inout) :: f_vbar(LBi:,LBj:)
322# endif
323# ifdef LCZ_FINAL
324 real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
325 real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
326# endif
327# endif
328 real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:)
329# ifdef HESSIAN_FSV
330 real(r8), intent(inout) :: f_zeta(LBi:,LBj:)
331# endif
332# ifdef LCZ_FINAL
333 real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:)
334# endif
335# else
336# ifdef MASKING
337 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
338 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
339 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
340# endif
341# ifdef ADJUST_BOUNDARY
342# ifdef SOLVE3D
343 real(r8), intent(inout) :: tl_t_obc(LBij:UBij,N(ng),4, &
344 & Nbrec(ng),2,NT(ng))
345 real(r8), intent(inout) :: tl_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
346 real(r8), intent(inout) :: tl_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
347# ifdef LCZ_FINAL
348 real(r8), intent(inout) :: ad_t_obc(LBij:UBij,N(ng),4, &
349 & Nbrec(ng),2,NT(ng))
350 real(r8), intent(inout) :: ad_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
351 real(r8), intent(inout) :: ad_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
352# endif
353# endif
354 real(r8), intent(inout) :: tl_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
355 real(r8), intent(inout) :: tl_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
356 real(r8), intent(inout) :: tl_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
357# ifdef LCZ_FINAL
358 real(r8), intent(inout) :: ad_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
359 real(r8), intent(inout) :: ad_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
360 real(r8), intent(inout) :: ad_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
361# endif
362# endif
363# ifdef ADJUST_WSTRESS
364 real(r8), intent(inout) :: tl_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
365 real(r8), intent(inout) :: tl_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
366# ifdef LCZ_FINAL
367 real(r8), intent(inout) :: ad_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
368 real(r8), intent(inout) :: ad_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
369# endif
370# endif
371# if defined ADJUST_STFLUX && defined SOLVE3D
372 real(r8), intent(inout) :: tl_tflux(LBi:UBi,LBj:UBj, &
373 & Nfrec(ng),2,NT(ng))
374# ifdef LCZ_FINAL
375 real(r8), intent(inout) :: ad_tflux(LBi:UBi,LBj:UBj, &
376 & Nfrec(ng),2,NT(ng))
377# endif
378# endif
379# ifdef SOLVE3D
380 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
381 real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
382 real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
383# if defined WEAK_CONSTRAINT && defined TIME_CONV
384 real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:)
385 real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:)
386# endif
387# ifdef HESSIAN_FSV
388 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
389 real(r8), intent(inout) :: f_t(LBi:UBi,LBj:UBj,N(ng),NT(ng))
390 real(r8), intent(inout) :: f_u(LBi:UBi,LBj:UBj,N(ng))
391 real(r8), intent(inout) :: f_v(LBi:UBi,LBj:UBj,N(ng))
392 real(r8), intent(inout) :: f_ubar(LBi:UBi,LBj:UBj)
393 real(r8), intent(inout) :: f_vbar(LBi:UBi,LBj:UBj)
394# endif
395# ifdef LCZ_FINAL
396 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
397 real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
398 real(r8), intent(inout) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
399# if defined WEAK_CONSTRAINT && defined TIME_CONV
400 real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:)
401 real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:)
402# endif
403# endif
404# else
405 real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:)
406 real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:)
407# ifdef HESSIAN_FSV
408 real(r8), intent(inout) :: f_ubar(LBi:UBi,LBj:UBj)
409 real(r8), intent(inout) :: f_vbar(LBi:UBi,LBj:UBj)
410# endif
411# ifdef LCZ_FINAL
412 real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:)
413 real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:)
414# endif
415# endif
416 real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,:)
417# ifdef HESSIAN_FSV
418 real(r8), intent(inout) :: f_zeta(LBi:UBi,LBj:UBj)
419# endif
420# ifdef LCZ_FINAL
421 real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:)
422# endif
423# endif
424
425
426
427 integer :: Lwrk, i, j, lstr, outLoop, rec, info
428
429 integer :: ndefLCZ = 1
430# ifdef LCZ_FINAL
431 integer :: ndefLZE = 1
432 integer :: ncidsav
433# endif
434# ifdef SOLVE3D
435 integer :: itrc, k
436# ifdef HESSIAN_FSV
437 real(r8) :: cff1, cff2
438 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
439 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC
440# endif
441# endif
442 real(r8) :: fac, fac1, fac2
443 real(r8) :: zbeta
444
445 real(r8), dimension(0:NstateVar(ng)) :: dot
446 real(r8), dimension(Ninner) :: DotProd
447 real(r8), dimension(Ninner) :: bvector
448 real(r8), dimension(Ninner) :: work
449# ifdef LCZ_FINAL
450 real(r8), dimension(Ninner,Ninner) :: GStemp
451 real(r8), dimension(Ninner,Ninner) :: GSsub
452 real(r8), dimension(Ninner) :: work1
453# endif
454
455# ifdef LCZ_FINAL
456 real :: sum
457 logical, save :: first = .true.
458 logical, save :: first1 = .true.
459# endif
460
461 character (len=256) :: ncname
462
463 character (len=*), parameter :: MyFile = &
464 & __FILE__//", tl_inner2state_tile"
465
466# include "set_bounds.h"
467
468 calledfrom=myfile
469 sourcefile=myfile
470
471
472
473
474
475
476 outloop=nouter
477 IF (master) THEN
478 DO i=1,ninner
479 zlanczos_diag(i)=cg_delta(i,outloop)
480 END DO
481 DO i=2,ninner
482 zlanczos_offdiag(i-1)=cg_beta(i,outloop)
483 END DO
484
485 CALL dpttrf (ninner, zlanczos_diag, zlanczos_offdiag, info)
486
487
488
489 DO i=1,ninner
490 zlanczos_diag(i)=sqrt(zlanczos_diag(i))
491 END DO
492 END IF
493# ifdef DISTRIBUTE
494 CALL mp_bcasti (ng, itlm, info)
495 CALL mp_bcastf (ng, itlm, zlanczos_diag)
496 CALL mp_bcastf (ng, itlm, zlanczos_offdiag)
497# endif
498 IF (info.ne.0) THEN
499 IF (master) WRITE (stdout,*) ' Error in DPTTRF: info = ', info
500 exit_flag=8
501 RETURN
502 END IF
503
504
505
506
507
508 DO i=1,ninner
509 work(i)=state(i)/zlanczos_diag(i)
510 END DO
511 DO i=ninner-1,1,-1
512 work(i)=work(i)-zlanczos_offdiag(i)*work(i+1)
513 END DO
514
515# ifdef LCZ_FINAL
516
517
518
519
520
521
522
523
524
525 SELECT CASE (lcz(ng)%IOtype)
526 CASE (io_nf90)
527 CALL netcdf_get_ivar (ng, itlm, trim(lcz(ng)%name), &
528 & 'ndefADJ', ndeflcz)
529
530# if defined PIO_LIB && defined DISTRIBUTE
531 CASE (io_pio)
532 CALL pio_netcdf_get_ivar (ng, itlm, trim(lcz(ng)%name), &
533 & 'ndefADJ', ndeflcz)
534# endif
535 END SELECT
536 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
537
538 IF (first) THEN
539 first=.false.
540
541
542
543
544
545
546
547
548 IF (ndeflcz.gt.0) THEN
549 lstr=len_trim(lcz(ng)%name)
550 WRITE (ncname,10) lcz(ng)%name(1:lstr-8), nouter
551 ELSE
552 ncname=lcz(ng)%name
553 END IF
554 lwrk=2
555
556
557
558 DO j=1,ninner
559 DO i=1,ninner
560 gsmatrix(i,j)=0.0_r8
561 END DO
562 END do3
563 DO j=1,ninner
564 gsmatrix(j,j)=1.0_r8
565 END DO
566
567 DO inner=1,ninner
568
569
570
571 DO j=1,ninner
572 DO i=1,ninner
573 gstemp(i,j)=0.0_r8
574 END DO
575 END DO
576 DO j=1,ninner
577 gstemp(j,j)=1.0_r8
578 END DO
579
580 ncname=lcz(ng)%name
581 CALL state_read (ng, tile, itlm, lcz(ng)%IOtype, &
582 & lbi, ubi, lbj, ubj, lbij, ubij, &
583 & lwrk, inner, &
584 & ndeflcz, lcz(ng)%ncid, &
585# if defined PIO_LIB && defined DISTRIBUTE
586 & lcz(ng)%pioFile, &
587# endif
588 & trim(ncname), &
589# ifdef MASKING
590 & rmask, umask, vmask, &
591# endif
592# ifdef ADJUST_BOUNDARY
593# ifdef SOLVE3D
594 & tl_t_obc, tl_u_obc, tl_v_obc, &
595# endif
596 & tl_ubar_obc, tl_vbar_obc, &
597 & tl_zeta_obc, &
598# endif
599# ifdef ADJUST_WSTRESS
600 & tl_ustr, tl_vstr, &
601# endif
602# if defined ADJUST_STFLUX && defined SOLVE3D
603 & tl_tflux, &
604# endif
605# ifdef SOLVE3D
606 & tl_t, tl_u, tl_v, &
607# else
608 & tl_ubar, tl_vbar, &
609# endif
610 & tl_zeta)
611
612
613
614
615
616
617
618
619 CALL state_copy (ng, tile, &
620 & lbi, ubi, lbj, ubj, lbij, ubij, &
621 & lwrk, lini, &
622# ifdef ADJUST_BOUNDARY
623# ifdef SOLVE3D
624 & tl_t_obc, tl_t_obc, &
625 & tl_u_obc, tl_u_obc, &
626 & tl_v_obc, tl_v_obc, &
627# endif
628 & tl_ubar_obc, tl_ubar_obc, &
629 & tl_vbar_obc, tl_vbar_obc, &
630 & tl_zeta_obc, tl_zeta_obc, &
631# endif
632# ifdef ADJUST_WSTRESS
633 & tl_ustr, tl_ustr, &
634 & tl_vstr, tl_vstr, &
635# endif
636# ifdef SOLVE3D
637# ifdef ADJUST_STFLUX
638 & tl_tflux, tl_tflux, &
639# endif
640 & tl_t, tl_t, &
641 & tl_u, tl_u, &
642 & tl_v, tl_v, &
643# if defined WEAK_CONSTRAINT && defined TIME_CONV
644 & tl_ubar, tl_ubar, &
645 & tl_vbar, tl_vbar, &
646# endif
647# else
648 & tl_ubar, tl_ubar, &
649 & tl_vbar, tl_vbar, &
650# endif
651 & tl_zeta, tl_zeta)
652
653
654
655 DO rec=1,inner-1
656
657
658
659
660 ncname=lze(ng)%name
661 CALL state_read (ng, tile, itlm, lze(ng)%IOtype, &
662 & lbi, ubi, lbj, ubj, lbij, ubij, &
663 & lwrk, rec, &
664 & ndeflze, lze(ng)%ncid, &
665# if defined PIO_LIB && defined DISTRIBUTE
666 & lze(ng)%pioFile, &
667# endif
668 & trim(ncname), &
669# ifdef MASKING
670 & rmask, umask, vmask, &
671# endif
672# ifdef ADJUST_BOUNDARY
673# ifdef SOLVE3D
674 & ad_t_obc, ad_u_obc, ad_v_obc, &
675# endif
676 & ad_ubar_obc, ad_vbar_obc, &
677 & ad_zeta_obc, &
678# endif
679# ifdef ADJUST_WSTRESS
680 & ad_ustr, ad_vstr, &
681# endif
682# if defined ADJUST_STFLUX && defined SOLVE3D
683 & ad_tflux, &
684# endif
685# ifdef SOLVE3D
686 & ad_t, ad_u, ad_v, &
687# else
688 & ad_ubar, ad_vbar, &
689# endif
690 & ad_zeta)
691 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
692
693
694
695 CALL state_dotprod (ng, tile, itlm, &
696 & lbi, ubi, lbj, ubj, lbij, ubij, &
697 & nstatevar(ng), dot(0:), &
698# ifdef MASKING
699 & rmask, umask, vmask, &
700# endif
701# ifdef ADJUST_BOUNDARY
702# ifdef SOLVE3D
703 & tl_t_obc(:,:,:,:,lini,:), &
704 & ad_t_obc(:,:,:,:,lwrk,:), &
705 & tl_u_obc(:,:,:,:,lini), &
706 & ad_u_obc(:,:,:,:,lwrk), &
707 & tl_v_obc(:,:,:,:,lini), &
708 & ad_v_obc(:,:,:,:,lwrk), &
709# endif
710 & tl_ubar_obc(:,:,:,lini), &
711 & ad_ubar_obc(:,:,:,lwrk), &
712 & tl_vbar_obc(:,:,:,lini), &
713 & ad_vbar_obc(:,:,:,lwrk), &
714 & tl_zeta_obc(:,:,:,lini), &
715 & ad_zeta_obc(:,:,:,lwrk), &
716# endif
717# ifdef ADJUST_WSTRESS
718 & tl_ustr(:,:,:,lini), ad_ustr(:,:,:,lwrk), &
719 & tl_vstr(:,:,:,lini), ad_vstr(:,:,:,lwrk), &
720# endif
721# ifdef SOLVE3D
722# ifdef ADJUST_STFLUX
723 & tl_tflux(:,:,:,lini,:), &
724 & ad_tflux(:,:,:,lwrk,:), &
725# endif
726 & tl_t(:,:,:,lini,:), ad_t(:,:,:,lwrk,:), &
727 & tl_u(:,:,:,lini), ad_u(:,:,:,lwrk), &
728 & tl_v(:,:,:,lini), ad_v(:,:,:,lwrk), &
729# else
730 & tl_ubar(:,:,lini), ad_ubar(:,:,lwrk), &
731 & tl_vbar(:,:,lini), ad_vbar(:,:,lwrk), &
732# endif
733 & tl_zeta(:,:,lini), ad_zeta(:,:,lwrk))
734
735 gstemp(rec,inner)=-dot(0)
736
737
738
739 fac1=1.0_r8
740 fac2=-dot(0)
741
742
743
744
745
746
747 CALL state_addition (ng, tile, &
748 & lbi, ubi, lbj, ubj, lbij, ubij, &
749 & lini, lwrk, lini, fac1, fac2, &
750# ifdef MASKING
751 & rmask, umask, vmask, &
752# endif
753# ifdef ADJUST_BOUNDARY
754# ifdef SOLVE3D
755 & tl_t_obc, ad_t_obc, &
756 & tl_u_obc, ad_u_obc, &
757 & tl_v_obc, ad_v_obc, &
758# endif
759 & tl_ubar_obc, ad_ubar_obc, &
760 & tl_vbar_obc, ad_vbar_obc, &
761 & tl_zeta_obc, ad_zeta_obc, &
762# endif
763# ifdef ADJUST_WSTRESS
764 & tl_ustr, ad_ustr, &
765 & tl_vstr, ad_vstr, &
766# endif
767# ifdef SOLVE3D
768# ifdef ADJUST_STFLUX
769 & tl_tflux, ad_tflux, &
770# endif
771 & tl_t, ad_t, &
772 & tl_u, ad_u, &
773 & tl_v, ad_v, &
774# if defined WEAK_CONSTRAINT && defined TIME_CONV
775 & tl_ubar, ad_ubar, &
776 & tl_vbar, ad_vbar, &
777# endif
778# else
779 & tl_ubar, ad_ubar, &
780 & tl_vbar, ad_vbar, &
781# endif
782 & tl_zeta, ad_zeta)
783 END DO
784
785
786
787 CALL state_dotprod (ng, tile, itlm, &
788 & lbi, ubi, lbj, ubj, lbij, ubij, &
789 & nstatevar(ng), dot(0:), &
790# ifdef MASKING
791 & rmask, umask, vmask, &
792# endif
793# ifdef ADJUST_BOUNDARY
794# ifdef SOLVE3D
795 & tl_t_obc(:,:,:,:,lini,:), &
796 & tl_t_obc(:,:,:,:,lini,:), &
797 & tl_u_obc(:,:,:,:,lini), &
798 & tl_u_obc(:,:,:,:,lini), &
799 & tl_v_obc(:,:,:,:,lini), &
800 & tl_v_obc(:,:,:,:,lini), &
801# endif
802 & tl_ubar_obc(:,:,:,lini), &
803 & tl_ubar_obc(:,:,:,lini), &
804 & tl_vbar_obc(:,:,:,lini), &
805 & tl_vbar_obc(:,:,:,lini), &
806 & tl_zeta_obc(:,:,:,lini), &
807 & tl_zeta_obc(:,:,:,lini), &
808# endif
809# ifdef ADJUST_WSTRESS
810 & tl_ustr(:,:,:,lini), tl_ustr(:,:,:,lini), &
811 & tl_vstr(:,:,:,lini), tl_vstr(:,:,:,lini), &
812# endif
813# ifdef SOLVE3D
814# ifdef ADJUST_STFLUX
815 & tl_tflux(:,:,:,lini,:), &
816 & tl_tflux(:,:,:,lini,:), &
817# endif
818 & tl_t(:,:,:,lini,:), tl_t(:,:,:,lini,:), &
819 & tl_u(:,:,:,lini), tl_u(:,:,:,lini), &
820 & tl_v(:,:,:,lini), tl_v(:,:,:,lini), &
821# else
822 & tl_ubar(:,:,lini), tl_ubar(:,:,lini), &
823 & tl_vbar(:,:,lini), tl_vbar(:,:,lini), &
824# endif
825 & tl_zeta(:,:,lini), tl_zeta(:,:,lini))
826
827 DO i=1,inner
828 gstemp(i,inner)=gstemp(i,inner)/sqrt(dot(0))
829 END DO
830
831
832
833
834
835 fac=1.0_r8/sqrt(dot(0))
836
837 CALL state_scale (ng, tile, &
838 & lbi, ubi, lbj, ubj, lbij, ubij, &
839 & lini, lini, fac, &
840# ifdef MASKING
841 & rmask, umask, vmask, &
842# endif
843# ifdef ADJUST_BOUNDARY
844# ifdef SOLVE3D
845 & tl_t_obc, tl_u_obc, tl_v_obc, &
846# endif
847 & tl_ubar_obc, tl_vbar_obc, &
848 & tl_zeta_obc, &
849# endif
850# ifdef ADJUST_WSTRESS
851 & tl_ustr, tl_vstr, &
852# endif
853# ifdef SOLVE3D
854# ifdef ADJUST_STFLUX
855 & tl_tflux, &
856# endif
857 & tl_t, tl_u, tl_v, &
858# else
859 & tl_ubar, tl_vbar, &
860# endif
861 & tl_zeta)
862
863
864
865
866
867
868
869
870
871 CALL state_copy (ng, tile, &
872 & lbi, ubi, lbj, ubj, lbij, ubij, &
873 & lini, lini, &
874# ifdef ADJUST_BOUNDARY
875# ifdef SOLVE3D
876 & ad_t_obc, tl_t_obc, &
877 & ad_u_obc, tl_u_obc, &
878 & ad_v_obc, tl_v_obc, &
879# endif
880 & ad_ubar_obc, tl_ubar_obc, &
881 & ad_vbar_obc, tl_vbar_obc, &
882 & ad_zeta_obc, tl_zeta_obc, &
883# endif
884# ifdef ADJUST_WSTRESS
885 & ad_ustr, tl_ustr, &
886 & ad_vstr, tl_vstr, &
887# endif
888# ifdef SOLVE3D
889# ifdef ADJUST_STFLUX
890 & ad_tflux, tl_tflux, &
891# endif
892 & ad_t, tl_t, &
893 & ad_u, tl_u, &
894 & ad_v, tl_v, &
895# if defined WEAK_CONSTRAINT && defined TIME_CONV
896 & ad_ubar, tl_ubar, &
897 & ad_vbar, tl_vbar, &
898# endif
899# else
900 & ad_ubar, tl_ubar, &
901 & ad_vbar, tl_vbar, &
902# endif
903 & ad_zeta, tl_zeta)
904
905
906
907
908 IF (inner.gt.1) THEN
909 lze(ng)%ncid=ncidsav
910 ELSE
911 ncidsav=lze(ng)%ncid
912 END IF
913 lze(ng)%Rindex=inner-1
914 CALL wrt_evolved (ng, lini, lini)
915 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
916
917
918
919
920 DO j=1,ninner
921 DO i=1,ninner
922 sum=0.0_r8
923 DO k=1,ninner
924 sum=sum+gsmatrix(i,k)*gstemp(k,j)
925 END DO
926 gssub(i,j)=sum
927 END DO
928 END DO
929 DO j=1,ninner
930 DO i=1,ninner
931 gsmatrix(i,j)=gssub(i,j)
932 END DO
933 END DO
934
935 END DO
936
937
938
939
940 IF (master) THEN
941 WRITE (stdout,*) ' Test of orthonormalization'
942 END IF
943
944 DO inner=1,ninner
945
946 ncname=lze(ng)%name
947 CALL state_read (ng, tile, itlm, lze(ng)%IOtype, &
948 & lbi, ubi, lbj, ubj, lbij, ubij, &
949 & lwrk, inner, &
950 & ndeflze, lze(ng)%ncid, &
951# if defined PIO_LIB && defined DISTRIBUTE
952 & lze(ng)%pioFile, &
953# endif
954 & trim(ncname), &
955# ifdef MASKING
956 & rmask, umask, vmask, &
957# endif
958# ifdef ADJUST_BOUNDARY
959# ifdef SOLVE3D
960 & tl_t_obc, tl_u_obc, tl_v_obc, &
961# endif
962 & tl_ubar_obc, tl_vbar_obc, &
963 & tl_zeta_obc, &
964# endif
965# ifdef ADJUST_WSTRESS
966 & tl_ustr, tl_vstr, &
967# endif
968# if defined ADJUST_STFLUX && defined SOLVE3D
969 & tl_tflux, &
970# endif
971# ifdef SOLVE3D
972 & tl_t, tl_u, tl_v, &
973# else
974 & tl_ubar, tl_vbar, &
975# endif
976 & tl_zeta)
977 DO rec=1,ninner
978
979 ncname=lze(ng)%name
980 CALL state_read (ng, tile, itlm, lze(ng)%IOtype, &
981 & lbi, ubi, lbj, ubj, lbij, ubij, &
982 & lwrk, rec, &
983 & ndeflze, lze(ng)%ncid, &
984# if defined PIO_LIB && defined DISTRIBUTE
985 & lze(ng)%pioFile, &
986# endif
987 & trim(ncname), &
988# ifdef MASKING
989 & rmask, umask, vmask, &
990# endif
991# ifdef ADJUST_BOUNDARY
992# ifdef SOLVE3D
993 & ad_t_obc, ad_u_obc, ad_v_obc, &
994# endif
995 & ad_ubar_obc, ad_vbar_obc, &
996 & ad_zeta_obc, &
997# endif
998# ifdef ADJUST_WSTRESS
999 & ad_ustr, ad_vstr, &
1000# endif
1001# if defined ADJUST_STFLUX && defined SOLVE3D
1002 & ad_tflux, &
1003# endif
1004# ifdef SOLVE3D
1005 & ad_t, ad_u, ad_v, &
1006# else
1007 & ad_ubar, ad_vbar, &
1008# endif
1009 & ad_zeta)
1010 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1011
1012
1013
1014 CALL state_dotprod (ng, tile, itlm, &
1015 & lbi, ubi, lbj, ubj, lbij, ubij, &
1016 & nstatevar(ng), dot(0:), &
1017# ifdef MASKING
1018 & rmask, umask, vmask, &
1019# endif
1020# ifdef ADJUST_BOUNDARY
1021# ifdef SOLVE3D
1022 & tl_t_obc(:,:,:,:,lwrk,:), &
1023 & ad_t_obc(:,:,:,:,lwrk,:), &
1024 & tl_u_obc(:,:,:,:,lwrk), &
1025 & ad_u_obc(:,:,:,:,lwrk), &
1026 & tl_v_obc(:,:,:,:,lwrk), &
1027 & ad_v_obc(:,:,:,:,lwrk), &
1028# endif
1029 & tl_ubar_obc(:,:,:,lwrk), &
1030 & ad_ubar_obc(:,:,:,lwrk), &
1031 & tl_vbar_obc(:,:,:,lwrk), &
1032 & ad_vbar_obc(:,:,:,lwrk), &
1033 & tl_zeta_obc(:,:,:,lwrk), &
1034 & ad_zeta_obc(:,:,:,lwrk), &
1035# endif
1036# ifdef ADJUST_WSTRESS
1037 & tl_ustr(:,:,:,lwrk), ad_ustr(:,:,:,lwrk), &
1038 & tl_vstr(:,:,:,lwrk), ad_vstr(:,:,:,lwrk), &
1039# endif
1040# ifdef SOLVE3D
1041# ifdef ADJUST_STFLUX
1042 & tl_tflux(:,:,:,lwrk,:), &
1043 & ad_tflux(:,:,:,lwrk,:), &
1044# endif
1045 & tl_t(:,:,:,lwrk,:), ad_t(:,:,:,lwrk,:), &
1046 & tl_u(:,:,:,lwrk), ad_u(:,:,:,lwrk), &
1047 & tl_v(:,:,:,lwrk), ad_v(:,:,:,lwrk), &
1048# else
1049 & tl_ubar(:,:,lwrk), ad_ubar(:,:,lwrk), &
1050 & tl_vbar(:,:,lwrk), ad_vbar(:,:,lwrk), &
1051# endif
1052 & tl_zeta(:,:,lwrk), ad_zeta(:,:,lwrk))
1053
1054 IF (master) THEN
1055 WRITE (stdout,*) 'inner = ', inner, ' rec = ', rec, &
1056 & ' dot-product = ', dot(0)
1057 END IF
1058
1059 END DO
1060 END DO
1061
1062
1063
1064 IF (master) THEN
1065 DO j=1,ninner
1066 DO i=1,ninner
1067 gsmatinv(i,j)=gsmatrix(i,j)
1068 END DO
1069 END DO
1070 CALL dtrtri ('U','N',ninner,gsmatinv,ninner,info)
1071 END IF
1072# ifdef DISTRIBUTE
1073 CALL mp_bcasti (ng, itlm, info)
1074 CALL mp_bcastf (ng, itlm, gsmatrix)
1075 CALL mp_bcastf (ng, itlm, gsmatinv)
1076# endif
1077 IF (info.ne.0) THEN
1078 IF (master) WRITE (stdout,*) ' Error in DPTTRF: info = ', info
1079 exit_flag=8
1080 RETURN
1081 END IF
1082
1083 END IF
1084
1085
1086
1087 DO i=1,ninner
1088 sum=0.0_r8
1089 DO j=1,ninner
1090 sum=sum+gsmatinv(i,j)*work(j)
1091 END DO
1092 work1(i)=sum
1093 END DO
1094 DO i=1,ninner
1095 work(i)=work1(i)
1096 END DO
1097
1098# endif
1099
1100
1101
1102
1103
1104
1105# ifdef LCZ_FINAL
1106
1107
1108 SELECT CASE (lze(ng)%IOtype)
1109 CASE (io_nf90)
1110 CALL netcdf_get_ivar (ng, itlm, trim(lze(ng)%name), &
1111 & 'ndefADJ', ndeflze)
1112
1113# if defined PIO_LIB && defined DISTRIBUTE
1114 CASE (io_pio)
1115 CALL pio_netcdf_get_ivar (ng, itlm, trim(lze(ng)%name), &
1116 & 'ndefADJ', ndeflze)
1117# endif
1118 END SELECT
1119 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1120# else
1121
1122
1123 SELECT CASE (lcz(ng)%IOtype)
1124 CASE (io_nf90)
1125 CALL netcdf_get_ivar (ng, itlm, trim(lcz(ng)%name), &
1126 & 'ndefADJ', ndeflcz)
1127
1128# if defined PIO_LIB && defined DISTRIBUTE
1129 CASE (io_pio)
1130 CALL pio_netcdf_get_ivar (ng, itlm, trim(lcz(ng)%name), &
1131 & 'ndefADJ', ndeflcz)
1132# endif
1133 END SELECT
1134 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1135# endif
1136
1137# ifdef LCZ_FINAL
1138
1139
1140
1141 IF (ndeflze.gt.0) THEN
1142 lstr=len_trim(lze(ng)%name)
1143 WRITE (ncname,10) lze(ng)%name(1:lstr-8), nouter
1144 10 FORMAT (a,'_',i4.4,'.nc')
1145 ELSE
1146 ncname=lze(ng)%name
1147 END IF
1148# else
1149
1150
1151
1152
1153
1154
1155
1156 IF (ndeflcz.gt.0) THEN
1157 lstr=len_trim(lcz(ng)%name)
1158 WRITE (ncname,10) lcz(ng)%name(1:lstr-8), nouter
1159 10 FORMAT (a,'_',i4.4,'.nc')
1160 ELSE
1161 ncname=lcz(ng)%name
1162 END IF
1163# endif
1164
1165 lwrk=2
1166# ifdef LCZ_FINAL
1167
1168
1169
1170
1171 IF (first1) THEN
1172 first1=.false.
1173 lze(ng)%ncid=ncidsav
1174 END IF
1175
1176# endif
1177 DO inner=1,ninner
1178
1179
1180
1181
1182
1183 CALL state_read (ng, tile, itlm, &
1184# ifdef LCZ_FINAL
1185 & lze(ng)%IOtype, &
1186# else
1187 & lcz(ng)%IOtype, &
1188# endif
1189 & lbi, ubi, lbj, ubj, lbij, ubij, &
1190 & lwrk, inner, &
1191# ifdef LCZ_FINAL
1192 & ndeflze, lze(ng)%ncid, &
1193# if defined PIO_LIB && defined DISTRIBUTE
1194 & lze(ng)%pioFile, &
1195# endif
1196# else
1197 & ndeflcz, lcz(ng)%ncid, &
1198# if defined PIO_LIB && defined DISTRIBUTE
1199 & lcz(ng)%pioFile, &
1200# endif
1201# endif
1202 & trim(ncname), &
1203# ifdef MASKING
1204 & rmask, umask, vmask, &
1205# endif
1206# ifdef ADJUST_BOUNDARY
1207# ifdef SOLVE3D
1208 & tl_t_obc, tl_u_obc, tl_v_obc, &
1209# endif
1210 & tl_ubar_obc, tl_vbar_obc, &
1211 & tl_zeta_obc, &
1212# endif
1213# ifdef ADJUST_WSTRESS
1214 & tl_ustr, tl_vstr, &
1215# endif
1216# if defined ADJUST_STFLUX && defined SOLVE3D
1217 & tl_tflux, &
1218# endif
1219# ifdef SOLVE3D
1220 & tl_t, tl_u, tl_v, &
1221# else
1222 & tl_ubar, tl_vbar, &
1223# endif
1224 & tl_zeta)
1225 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1226
1227 IF (inner.eq.1) THEN
1228 fac=0.0_r8
1229 CALL state_initialize (ng, tile, &
1230 & lbi, ubi, lbj, ubj, lbij, ubij, &
1231 & lini, fac, &
1232# ifdef MASKING
1233 & rmask, umask, vmask, &
1234# endif
1235# ifdef ADJUST_BOUNDARY
1236# ifdef SOLVE3D
1237 & tl_t_obc, tl_u_obc, tl_v_obc, &
1238# endif
1239 & tl_ubar_obc, tl_vbar_obc, &
1240 & tl_zeta_obc, &
1241# endif
1242# ifdef ADJUST_WSTRESS
1243 & tl_ustr, tl_vstr, &
1244# endif
1245# if defined ADJUST_STFLUX && defined SOLVE3D
1246 & tl_tflux, &
1247# endif
1248# ifdef SOLVE3D
1249 & tl_t, tl_u, tl_v, &
1250# else
1251 & tl_ubar, tl_vbar, &
1252# endif
1253 & tl_zeta)
1254 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1255 END IF
1256
1257
1258
1259
1260
1261 fac1=1.0_r8
1262 fac2=work(inner)
1263 CALL state_addition (ng, tile, &
1264 & lbi, ubi, lbj, ubj, lbij, ubij, &
1265 & lini, lwrk, lini, &
1266 & fac1, fac2, &
1267# ifdef MASKING
1268 & rmask, umask, vmask, &
1269# endif
1270# ifdef ADJUST_BOUNDARY
1271# ifdef SOLVE3D
1272 & tl_t_obc, tl_t_obc, &
1273 & tl_u_obc, tl_u_obc, &
1274 & tl_v_obc, tl_v_obc, &
1275# endif
1276 & tl_ubar_obc, tl_ubar_obc, &
1277 & tl_vbar_obc, tl_vbar_obc, &
1278 & tl_zeta_obc, tl_zeta_obc, &
1279# endif
1280# ifdef ADJUST_WSTRESS
1281 & tl_ustr, tl_ustr, &
1282 & tl_vstr, tl_vstr, &
1283# endif
1284# if defined ADJUST_STFLUX && defined SOLVE3D
1285 & tl_tflux, tl_tflux, &
1286# endif
1287# ifdef SOLVE3D
1288 & tl_t, tl_t, &
1289 & tl_u, tl_u, &
1290 & tl_v, tl_v, &
1291# if defined WEAK_CONSTRAINT && defined TIME_CONV
1292 & tl_ubar, tl_ubar, &
1293 & tl_vbar, tl_vbar, &
1294# endif
1295# else
1296 & tl_ubar, tl_ubar, &
1297 & tl_vbar, tl_vbar, &
1298# endif
1299 & tl_zeta, tl_zeta)
1300 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1301
1302 END DO
1303# ifdef HESSIAN_FSV
1304
1305
1306
1307
1308 DO j=jstrr,jendr
1309 DO i=istrr,iendr
1310 f_zeta(i,j)=tl_zeta(i,j,lini)
1311 END DO
1312 END DO
1313# ifndef SOLVE3D
1314
1315
1316
1317 DO j=jstrr,jendr
1318 DO i=istr,iendr
1319 f_ubar(i,j)=tl_ubar(i,j,lini)
1320 END DO
1321 END DO
1322
1323 DO j=jstr,jendr
1324 DO i=istrr,iendr
1325 f_vbar(i,j)=tl_vbar(i,j,lini)
1326 END DO
1327 END DO
1328
1329# else
1330
1331
1332
1333 DO k=1,n(ng)
1334 DO j=jstrr,jendr
1335 DO i=istr,iendr
1336 f_u(i,j,k)=tl_u(i,j,k,lini)
1337 END DO
1338 END DO
1339 DO j=jstr,jendr
1340 DO i=istrr,iendr
1341 f_v(i,j,k)=tl_v(i,j,k,lini)
1342 END DO
1343 END DO
1344 END DO
1345
1346
1347
1348 DO j=jstr,jend
1349 DO i=istru,iend
1350 dc(i,0)=0.0_r8
1351 cf(i,0)=0.0_r8
1352 END DO
1353 DO k=1,n(ng)
1354 DO i=istru,iend
1355 dc(i,k)=0.5_r8*(hz(i,j,k)+hz(i-1,j,k))
1356 dc(i,0)=dc(i,0)+dc(i,k)
1357 cf(i,0)=cf(i,0)+dc(i,k)*f_u(i,j,k)
1358 END DO
1359 END DO
1360 DO i=istru,iend
1361 cff1=1.0_r8/dc(i,0)
1362 cff2=cf(i,0)*cff1
1363# ifdef MASKING
1364 cff2=cff2*umask(i,j)
1365# endif
1366 f_ubar(i,j)=cff2
1367 END DO
1368 END DO
1369
1370
1371
1372 DO j=jstrv,jend
1373 IF (j.ge.jstrm) THEN
1374 DO i=istr,iend
1375 dc(i,0)=0.0_r8
1376 cf(i,0)=0.0_r8
1377 END DO
1378 DO k=1,n(ng)
1379 DO i=istr,iend
1380 dc(i,k)=0.5_r8*(hz(i,j,k)+hz(i,j-1,k))
1381 dc(i,0)=dc(i,0)+dc(i,k)
1382 cf(i,0)=cf(i,0)+dc(i,k)*f_v(i,j,k)
1383 END DO
1384 END DO
1385 DO i=istr,iend
1386 cff1=1.0_r8/dc(i,0)
1387 cff2=cf(i,0)*cff1
1388# ifdef MASKING
1389 cff2=cff2*vmask(i,j)
1390# endif
1391 f_vbar(i,j)=cff2
1392 END DO
1393 END IF
1394 END DO
1395
1396
1397
1398 DO itrc=1,nt(ng)
1399 DO k=1,n(ng)
1400 DO j=jstrr,jendr
1401 DO i=istrr,iendr
1402 f_t(i,j,k,itrc)=tl_t(i,j,k,lini,itrc)
1403 END DO
1404 END DO
1405 END DO
1406 END DO
1407# endif
1408 fac=0.0_r8
1409 CALL state_initialize (ng, tile, &
1410 & lbi, ubi, lbj, ubj, lbij, ubij, &
1411 & lini, fac, &
1412# ifdef MASKING
1413 & rmask, umask, vmask, &
1414# endif
1415# ifdef ADJUST_BOUNDARY
1416# ifdef SOLVE3D
1417 & tl_t_obc, tl_u_obc, tl_v_obc, &
1418# endif
1419 & tl_ubar_obc, tl_vbar_obc, &
1420 & tl_zeta_obc, &
1421# endif
1422# ifdef ADJUST_WSTRESS
1423 & tl_ustr, tl_vstr, &
1424# endif
1425# if defined ADJUST_STFLUX && defined SOLVE3D
1426 & tl_tflux, &
1427# endif
1428# ifdef SOLVE3D
1429 & tl_t, tl_u, tl_v, &
1430# else
1431 & tl_ubar, tl_vbar, &
1432# endif
1433 & tl_zeta)
1434
1435# endif
1436
1437 RETURN