source: trunk/LMDZ.GENERIC/libf/phystd/newtrelax.F90 @ 1145

Last change on this file since 1145 was 787, checked in by aslmd, 12 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

  1. ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
  2. if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

  • Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
  • RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
  • RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
  • RUN RCM1D : perfect match in normal mode.
  • RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

  • Additional arguments to some subroutines (ngrid and nq)
  • F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

  • F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
  • Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
  • Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
  • Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

  • Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

  • nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this
File size: 3.5 KB
Line 
1subroutine newtrelax(ngrid,mu0,sinlat,popsk,temp,pplay,pplev,dtrad,firstcall)
2       
3  implicit none
4
5#include "dimensions.h"
6#include "dimphys.h"
7#include "comcstfi.h"
8#include "callkeys.h"
9#include "netcdf.inc"
10
11!==================================================================
12!     
13!     Purpose
14!     -------
15!     Alternative Newtonian radiative transfer scheme.
16!     
17!     Authors
18!     -------
19!     R. Wordsworth (2010)
20!     
21!==================================================================
22 
23  integer ngrid
24 
25  ! Input
26  real mu0(ngrid)                    ! cosine of sun incident angle
27  real sinlat(ngrid)                 ! sine of latitude
28  real temp(ngrid,nlayermx)          ! temperature at each layer (K)
29  real pplay(ngrid,nlayermx)         ! pressure at each layer (Pa)
30  real pplev(ngrid,nlayermx+1)       ! pressure at each level (Pa)
31  real popsk(ngrid,nlayermx)         ! pot. T to T converter
32
33  ! Output
34  real dtrad(ngrid,nlayermx)
35
36  ! Internal
37  real Trelax_V, Trelax_H
38  real,allocatable,dimension(:,:),save :: Trelax
39
40  real T_trop ! relaxation temperature at tropopause (K)
41  real T_surf ! relaxation temperature at surface (K)
42  real dT_EP  ! Equator-Pole relaxation temperature difference (K)
43
44  real sig, f_sig, sig_trop
45  integer l,ig
46  logical firstcall
47
48
49  logical tidallocked
50  parameter (tidallocked = .true.)
51
52  ! Setup relaxation temperature 
53  if(firstcall) then
54
55     ALLOCATE(Trelax(ngrid,nlayermx))
56
57     print*,'-----------------------------------------------------'
58     print*,'| ATTENTION: You are using a Newtonian cooling scheme'
59     print*,'| for the radiative transfer. This means that ALL'
60     print*,'| other physics subroutines must be switched off.'
61     print*,'-----------------------------------------------------'
62
63     if(tidallocked)then
64        do ig=1,ngrid
65
66           T_surf = 126. + 239.*mu0(ig)
67           T_trop = 140. + 52.*mu0(ig)
68           do l=1,nlayermx
69
70              if(mu0(ig).le.0.0)then ! night side
71                 Trelax(ig,l)=0.0
72              else                   ! day side
73                 Trelax(ig,l) = T_surf*popsk(ig,l)
74                 if (Trelax(ig,l).lt.T_trop) Trelax(ig,l) = T_trop
75              endif
76
77           enddo
78        enddo
79
80     else
81
82        T_trop = 200.
83        T_surf = 288.
84        dT_EP  = 70.
85
86        sig_trop=(T_trop/T_surf)**(1./rcp)
87
88        do l=1,nlayermx
89           do ig=1,ngrid
90
91              ! vertically varying component
92              Trelax_V = T_surf*popsk(ig,l)
93              if (Trelax_V.lt.T_trop) Trelax_V = T_trop
94             
95              ! horizontally varying component
96              sig = pplay(ig,l)/pplev(ig,1)
97              if(sig.ge.sig_trop)then
98                 f_sig=sin((pi/2)*((sig-sig_trop)/(1-sig_trop)))
99              else
100                 f_sig=0.0
101              endif
102              Trelax_H = -f_sig*dT_EP*(sinlat(ig)**2 - 1./3.)
103             
104              Trelax(ig,l) = Trelax_V + Trelax_H           
105           
106           enddo
107        enddo
108
109     endif
110
111     firstcall=.false.
112 
113  endif
114
115  ! Calculate radiative forcing
116  do l=1,nlayermx
117     do ig=1,ngrid
118        dtrad(ig,l) = -(temp(ig,l) - Trelax(ig,l)) / tau_relax
119        if(temp(ig,l).gt.500.)then ! Trelax(ig,l))then
120           print*,'ig=',ig
121           print*,'l=',l
122           print*,'temp=',temp(ig,l)
123           print*,'Trelax=',Trelax(ig,l)
124        endif
125     enddo
126  enddo
127
128  call writediagfi(ngrid,'Tref','rad forc temp','K',3,Trelax)
129  !call writediagfi(ngrid,'ThetaZ','stellar zenith angle','deg',2,mu0)
130
131  return
132end subroutine newtrelax
Note: See TracBrowser for help on using the repository browser.